fftstokes.m

% geoid computation using planar FFT by Stokes' formula

load grav80res.grd; %        
% GRS80   EIGEN-CG03C     
% ( 5)  

load geoid80gm.grd; %       
%  EIGEN-CG03C ( 5)

load geoid80top.grd; %      
%  ( 5) 

dx=5.28; %   
dy=9.27; %   
nrow=49; %   
ncol=73; %   

%         . 
%     () cx    (/2)+1   cx
%    .       cx 
%     .  cx     
%  .
if mod(nrow,2)==0 %mod:  
    cx=nrow/2+1;
else
    cx=round(nrow/2); %    
end

%            
%   .     
%() cy         .
%  cx    cy       
% .
if mod(ncol,2)==0
    cy=ncol/2+1;
else
    cy=round(ncol/2);
end

%         
%  l_N=(r^(-1) (x,y))^2.
for i=nrow:-1:1 %   i  49-1 
for j=1:1:ncol  %   j  1-73
%   i,j   1/r^2 (i,j)     
%   kernel.   abs    
%    cy-j  cx-i.    
%  dx  dy .    
%     r^2 (i,j    
kernel(i,j)=1/((dx*abs(cy-j))^2+(dy*abs(cx-i))^2);
%           
%   .     r^2 (i,j)     
%      kernel        
%.
if  kernel(i,j)==inf
    kernel(i,j)=0;
end
end
end
%     l_N=r^(-1) (x,y) 
kernels=sqrt(kernel);
%   C=R/2_m
C=1000./(2*pi*979000.0);
%    _g,l_N    
% Fourier 
DG=fft2(grav80res);
kernel2=fft2(kernels);
%      
%   f1   -1
f1=-1;
for i=nrow:-1:1 %   i  49-1
 for j=1:1:ncol %   j  1-73
     %   (i,j)   ?  
     %    DG  kernel2.     
     %    1,  i+j ,  -1,  i+j 
     %.
        N(i,j)=DG(i,j).*kernel2(i,j)*f1^(i+j);
    end
end
%    DN    
% Fourier   . (     
%    ifft(F_1 (t) F_2 (t))  ifft() 
%  Fourier   F()   
%Fourier)
DN=ifft2(N);
%      DN    
%       
%real(DN)  C       Y.
NNN=real(DN).*C*dx*dy;
%         
%        
%       EIGEN-CG03C.   
for i=nrow:-1:1
 for j=1:1:ncol
        geoidfin(i,j)=geoid80gm(i,j)+geoid80top(i,j)+ NNN(i,j); %DG(i,j).*kernel2(i,j)*f1^(i+j);
 end
end
%        nfinal.dat  
%  5  . 
fid=fopen('nfinal.dat','w');
for i=1:1:nrow
 for j=1:1:ncol
    fprintf(fid,'%15.3f',geoidfin(i,j));
 end
 fprintf(fid,'\n');
end
fclose(fid);

%             
%      (1/12  )
f=60:(-1/12):56;
l=-120:(1/12):-114;

%    .    
% subplot        
% ,        .
%    xlabel,ylabel,title  
%      X   Y  
%  .        
%        
%  .
subplot (3,1,1)
contourf(l,f,NNN),title('RESIDUAL GEOID HEIGHTS')
xlabel('Longitude (deg)'),ylabel('Latitude (deg)'),colorbar
subplot (3,1,2)
contourf(l,f,geoid80gm),title('EIGEN-CG03C GEOID HEIGHTS')
xlabel('Longitude (deg)'),ylabel('Latitude (deg)'),colorbar
subplot(3,1,3)
contourf(l,f,geoid80top),title('REMAINING TOPOGRAPHY')
xlabel('Longitude (deg)'),ylabel('Latitude (deg)'),colorbar
figure
contourf(l,f,grav80res),title('RESIDUAL GRAVITY ANOMALIES')
xlabel('Longitude (deg)'),ylabel('Latitude (deg)'),colorbar
figure
contourf(l,f,geoidfin),title('FINAL GEOID HEIGHTS')
xlabel('Longitude (deg)'),ylabel('Latitude (deg)'),colorbar