function dgeo=dgeodesic(fx,fy)
% computes the geodesic distance between two power spectra

dt=.01; dgeo=0;
for t=0:dt:1-dt,
    fxy_t=fx.*(fy./fx).^t;
    fxy_tdt=fx.*(fy./fx).^(t+dt);
   dgeo=dgeo+sqrt(ag(fxy_t,fxy_tdt));
end