%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