Parker反演

%3DINVER.M: A MATLAB program to invert the gravity anomaly over a 3-D horizontal density interface by Parker-Oldenburg‘s algorithm
%David Gomez Ortiz and Bhrigu N P Agarwal
%
%Iterative process of Oldenburg（1974）
%using the Parker‘s （1973） equation
%Given a gravity map as input （bouin）in mGal
%as well as density contrast （contrast） in gr/cm3
%and the mean source depth （z0） in Km
%the program computes the topography of the interface （topoout） in Km 
%
%density contrast is positive （e. g. mantle-crust）
%mean depth reference is positive downwards（z0）
%It is also necessary to include map length in Km （longx and longy）
%considering a square area the number of rows or columns （numrows and numcolumns）
%and the convergence criterion （criterio） in Km
%
%The maximum number of iterations is 10
%
%The program computes also the gravity map （bouinv） due to the computed relief（forward modelling）
%The cut-off frequencies （SH and WH） for the filter in order to achieve
%convergence have to be also previously defined （as 1/wavelength in Km）
%

bouin=‘c:\gravityinput.dat‘;   %path and name of the input gravity file e.g. c:\gravityinput.dat
topoout=‘c:\topooutput.dat‘;  %path and name of the output topography file e.g. c:\topooutput.dat
bouinverted=‘c:\bouinv.dat‘;   %path and name of the output gravity file e.g. c:\bouinv.dat
numrows=102;                  %e. g.  112 data in Y direction
numcolumns=102;           % e. g. 112 data in X direction
longx=1212;              %e.g. 666 Km data length in X direction
longy=1212;            %e.g. 666 Km data length in Y direction
contrast=0.4;          %density contrast value in g/cm3 e. g. 0.4
criterio=0.02;		%convergence criterion in Km e. g. 0.02
z0=30;               %mean reference depth in Km e. g. 30 Km
WH=0.002;                %smaller cut-off frequency （1/wavelength in Km） e. g. 0.01 i.e. 1/100 Km
SH=0.0022;                %greater cut-off frequency （1/wavelength in Km） e. g. 0.012 i.e. 1/83.3 Km

truncation=0.1          %truncation window data length in % per one i.e 10%

fid=fopen（bouin）;
bou=fscanf（fid‘%f‘[numcolumns numrows]）;    %open the input gravity file read it and stores it in variable ‘bou‘
fclose（fid）

% Calculate mean value of the gravity map
fftbou=fft2（bou）;   %fft2 computes the 2-D FFT of a 2-D matrix （bou in this case）
meangravity=（fftbou（11）/（numrows*numcolumns））; %the first term of the 2-D fft array divided by the product of the number of rows and columns provides the mean gravity value of the map
  
% Demean data
disp（‘Data sets demeaned‘）  %displays the text ‘Data sets demeaned‘ on the screen
bou=bou-meangravity;  %input gravity map is demeaned

%A cosine Tukey window with a truncation of 10% default is applied
wrows = tukeywin（numrowstruncation）;  %this computes a 1-D cosine Tukey window of the same length as the original matrix input rows and with the truncation defined by the variable ‘truncation‘
wcolumns = tukeywin（numcolumns