function tercor=parker(c11,cmn,topo,ninterf,addm,drho) % tercor=PARKER(c11,cmn,topo,ninterf,addm,drho) % % Computes terrain/Bouguer corrections using Parker's method. % Bathymetry input is negative. The resulting correction must be ADDED to % the free-air anomaly. Everything in standard units (SI)! % So 'tercor' is in m/s2. To convert to milliGal, divide by 1e-5. % Take any nans' out, as in % tdint(isnan(tdint))=mean(tdint(~isnan(tdint))); % See Parker, GJRAS 1972 (31) 447-455 % Compare this to adding 2*pi*drho* the water depth to "replace" % the water with mantle density. % % INPUT: % % c11,cmn Coordinates of top left and bottom right pixel [degrees] % topo Topography/bathymetry matrix [m], % ninterf Number of interfaces % addm Depth to the interface(s) in [m>0] % drho Density contrast(s) [kg/m3] % % tercor=parker(c11,cmn,topo,2,[0 6000],[1700 600]); % % Last modified by fjsimons-at-alum.mit.edu, October 23rd, 2003 % Contributions from Mark Behn, Ban-Yuan Kuo and Jian Lin %---Define Input Parameters--- npower=5; % Power of Taylor Series Expansion ifold=1; % Folding flag (ifold=1 : folding, ifold=2 : no folding) %---Load Datafile--- [ny,nx]=size(topo); lenx=abs(cmn(1)-c11(1))*fralmanac('DegDis'); leny=abs(c11(2)-cmn(2))*fralmanac('DegDis'); dx=lenx/(nx-1); dy=leny/(ny-1); %---Normalize data to new reference level----- zmax=max(topo(:)); zmin=min(topo(:)); slev=mean(topo(:)); ddepth=slev-topo; %\$ How much deeper than mean depth % Gravitational constant in SI grav=2*pi*drho*fralmanac('GravCst'); if ifold==1 ddepth=flipflop(ddepth); lenx=2*lenx; leny=2*leny; end % Wave numbers rad/m K=knum2(size(ddepth),[leny lenx]); % Because fft is not fftshifted here not either kwn=fftshift(K); % Parker (1972) Eq. 4 csum=zeros(size(ddepth)); for ip=1:npower fprintf('%5s %12.4f\n','Power =',ip) data=ddepth.^ip; data=fft2(data); data=data.*(kwn.^(ip-1))/factorial(ip); csum=csum+data; end [k,l]=find(kwn==0); csum(k,l)=0; fprintf('%10s\n','Upward Continuation') tercor=zeros(ny,nx); for iref=1:ninterf % Depth over which to be upward continued, starting from mean elevation zlev=abs(slev)+addm(iref); % Upward continuation: the field must be attenuated data=csum*grav(iref).*exp(-zlev*kwn); fprintf('%10s\n','Inverse Transformation') data=ifft2(data); if mean(abs(imag(data(:))))>1e-14 warning(sprintf('Transformed data is not real by %s',... mean(abs(imag(data(:)))))) end ddepth=real(data(1:ny,1:nx)); tercor=tercor+ddepth; end