<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<HTML><HEAD>
<META http-equiv=Content-Type content="text/html; charset=iso-8859-1">
<META content="MSHTML 6.00.6000.16674" name=GENERATOR>
<STYLE></STYLE>
</HEAD>
<BODY bgColor=#ffffff>
<DIV><FONT face="Comic Sans MS">Hello</FONT></DIV>
<DIV><FONT face="Comic Sans MS"></FONT> </DIV>
<DIV><FONT face="Comic Sans MS">I have a Matlab file I would like to convert to 
sci (scilab) format but it fails to work. I am using scilab 4.1.2</FONT></DIV>
<DIV><FONT face="Comic Sans MS"></FONT> </DIV>
<DIV><FONT face="Comic Sans MS">I am not a Matlab user so cannot say why it is 
failing:</FONT></DIV>
<DIV><FONT face="Comic Sans MS"></FONT> </DIV>
<DIV><FONT face="Comic Sans MS">I get the following error:</FONT></DIV>
<DIV><FONT face="Comic Sans MS"></FONT> </DIV>
<DIV><FONT face="Comic Sans MS">Startup execution:<BR>  loading initial 
environment</FONT></DIV>
<DIV> </DIV>
<DIV><FONT face="Comic Sans MS">-->mfile2sci();<BR>  ****** Beginning of 
mfile2sci() session ******<BR>  File to convert: 
D:/3D_Inver/3dinver.m<BR>  Result file path: D:/Temp/<BR>  Recursive 
mode: OFF<BR>  Only double values used in M-file: NO<BR>  Verbose 
mode: 3<BR>  Generate formated code: NO<BR>  M-file 
reading...<BR>  M-file reading: Done<BR>  Syntax 
modification...<BR>  Syntax modification: Done<BR> !--error 
37<BR>incorrect function at line      1<BR>at 
line     190 of function mfile2sci called by 
:<BR>line    64 of function m2sci_gui called by 
:<BR>line    16 of function mfile2sci called by 
:<BR>mfile2sci();</FONT></DIV>
<DIV> </DIV><FONT face="Comic Sans MS">
<DIV><BR>--></DIV>
<DIV> </DIV>
<DIV>I have listed the file here in case someone can see what the issue is</DIV>
<DIV> </DIV>
<DIV>Thanks</DIV>
<DIV> </DIV>
<DIV>Lester</DIV>
<DIV> </DIV>
<DIV>%3DINVER.M: A MATLAB program to invert the gravity anomaly over a 3-D 
horizontal density interface by Parker-Oldenburg's algorithm<BR>%David Gomez 
Ortiz and Bhrigu N P Agarwal<BR>%<BR>%Iterative process of 
Oldenburg(1974)<BR>%using the Parker's (1973) equation<BR>%Given a gravity map 
as input (bouin)in mGal<BR>%as well as density contrast (contrast) in 
gr/cm3<BR>%and the mean source depth (z0) in Km<BR>%the program computes the 
topography of the interface (topoout) in Km <BR>%<BR>%density contrast is 
positive (e. g. mantle-crust)<BR>%mean depth reference is positive 
downwards(z0)<BR>%It is also necessary to include map length in Km (longx and 
longy)<BR>%considering a square area, the number of rows or columns (numrows and 
numcolumns)<BR>%and the convergence criterion (criterio) in Km<BR>%<BR>%The 
maximum number of iterations is 10<BR>%<BR>%The program computes also the 
gravity map (bouinv) due to the computed relief(forward modelling)<BR>%The 
cut-off frequencies (SH and WH) for the filter in order to 
achieve<BR>%convergence have to be also previously defined (as 1/wavelength in 
Km)<BR>%</DIV>
<DIV> </DIV>
<DIV>bouin='c:\gravityinput.dat';   %path and name of the input 
gravity file, e.g. c:\gravityinput.dat<BR>topoout='c:\topooutput.dat';  
%path and name of the output topography file, e.g. 
c:\topooutput.dat<BR>bouinverted='c:\bouinv.dat';   %path and name of 
the output gravity file, e.g. 
c:\bouinv.dat<BR>numrows=112;                  
%e. g.  112 data in Y 
direction<BR>numcolumns=112;           
% e. g. 112 data in X 
direction<BR>longx=666;              
%e.g. 666 Km data length in X 
direction<BR>longy=666;            
%e.g. 666 Km data length in Y 
direction<BR>contrast=0.4;          
%density contrast value in g/cm3, e. g. 
0.4<BR>criterio=0.02;  %convergence criterion in Km, e. g. 
0.02<BR>z0=30;               
%mean reference depth in Km, e. g. 30 
Km<BR>WH=0.01;                
%smaller cut-off frequency (1/wavelength in Km) e. g. 0.01, i.e. 1/100 
Km<BR>SH=0.012;                
%greater cut-off frequency (1/wavelength in Km) e. g. 0.012, i.e. 1/83.3 
Km</DIV>
<DIV> </DIV>
<DIV>truncation=0.1          
%truncation window data length in % per one, i.e, 10%</DIV>
<DIV> </DIV>
<DIV>fid=fopen(bouin);<BR>bou=fscanf(fid,'%f',[numcolumns 
numrows]);    %open the input gravity file, read it and stores it 
in variable 'bou'<BR>fclose(fid)</DIV>
<DIV> </DIV>
<DIV>% Calculate mean value of the gravity map<BR>fftbou=fft2(bou);   
%fft2 computes the 2-D FFT of a 2-D matrix (bou, in this 
case)<BR>meangravity=(fftbou(1,1)/(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<BR>  <BR>% Demean data<BR>disp('Data sets 
demeaned')  %displays the text 'Data sets demeaned' on the 
screen<BR>bou=bou-meangravity;  %input gravity map is demeaned</DIV>
<DIV> </DIV>
<DIV>%A cosine Tukey window with a truncation of 10% default is applied<BR>wrows 
= tukeywin(numrows,truncation);  %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'<BR>wcolumns = 
tukeywin(numcolumns,truncation);  %this computes a 1-D cosine Tukey window 
of the same length as the original matrix input columns and with the truncation 
defined by the variable 'truncation'</DIV>
<DIV> </DIV>
<DIV>w2 =wrows(:)* wcolumns(:).'; %this generates a 2-D cosine Tukey window 
multipliying the 1-D windows<BR>bou=bou.*w2';   %the original gravity 
input matrix, previously demeaned, is multiplied by the cosine window</DIV>
<DIV> </DIV>
<DIV>mapabou=bou';   %the original gravity input matrix after 
demeaning is transposed<BR>surf(mapabou)   %then it is drawn only for 
visualization<BR>title('Observed gravity anomaly map (mGal)') %this is the title 
of the graph</DIV>
<DIV> </DIV>
<DIV>fftbou=fft2(mapabou);  %the 2-D FFT of the gravity input matrix is 
computed after demeaning<BR>spectrum=abs(fftbou(1:numrows/2, 
1:numcolumns/2));  %this computes the amplitude 
spectrum<BR>figure<BR>surf(spectrum)   %and the spectrum is drawn only 
for visualization<BR>title('Amplitude spectrum of the Observed gravity map') 
%this is the title of the new graph</DIV>
<DIV> </DIV>
<DIV>%the matrix with the frequencies of every harmonic is computed<BR>for 
f=1:abs((numrows/2)+1);<BR>   for 
g=1:abs((numcolumns/2)+1);<BR>      frequency(f, 
g)=sqrt(((f-1)/longx)^2+((g-1)/longy)^2);<BR>   end;<BR>end;<BR>%the 
matrix of the negative frequencies is also 
computed<BR>frequency2=fliplr(frequency);<BR>frequency3=flipud(frequency);<BR>frequency4=fliplr(flipud(frequency));<BR>entero=round(numcolumns/2);<BR>if 
((numcolumns/2)- entero)==0<BR>   frequency2(:,1)=[];<BR>   
frequency3(1,:)=[];<BR>   frequency4(:,1)=[];<BR>   
frequency4(1,:)=[];<BR>   frequencytotal=[frequency 
frequency2;frequency3 frequency4];<BR>else<BR>   
frequencytotal=[frequency frequency2;frequency3 
frequency4];<BR>end<BR>frequencytotal(end,:)=[];<BR>frequencytotal(:,end)=[];</DIV>
<DIV> </DIV>
<DIV>frequencytotal=frequencytotal.*(2*pi);  %the frequency (1/wavelength) 
matrix is transformed to wavenumber (2*pi/wavelength) matrix</DIV>
<DIV> </DIV>
<DIV>%The iterative process starts here<BR>%The first term of the series, that 
is constant, is computed and stored in variable 
'constant'<BR>up=-(fftbou.*(exp((z0*100000)*(frequencytotal.*(1/100000)))));<BR>down=(2*pi*6.67*contrast);<BR>constant=up./down;<BR>%The 
high-cut filter is 
constructed<BR>filter=frequencytotal.*0;      %the 
filter matrix is set to zero<BR>frequencytotal=frequencytotal./(2*pi);  
<BR>for f=1:numrows;<BR>   for 
g=1:numcolumns;<BR>      if 
frequencytotal(f,g)<WH<BR>      filter(f,g)=1;  
<BR>elseif frequencytotal(f,g)<SH<BR>      
filter(f,g)=0.5.*(1+cos((((2*pi)*frequencytotal(f,g))-(2*pi*WH))/(2*(SH-WH))));<BR>else<BR>filter(f,g)=0;<BR>end<BR>   
end;<BR>end;<BR>constant=constant.*filter;  %the filter is applied to the 
first term of the series</DIV>
<DIV> </DIV>
<DIV>topoinverse=real(ifft2(constant));  %the real part of the inverse 2-D 
FFT of the first term provides the first topography approach stored in the 
variable 'topoinverse'<BR>frequencytotal=frequencytotal.*(2*pi);</DIV>
<DIV> </DIV>
<DIV>%It starts the computation of the second term of the series with a maximum 
of 10 iterations<BR>   
<BR>topo2=(((frequencytotal.^(1))./(prod(1:2))).*(fft2(topoinverse.^2)));   
%the function 'prod' computes the factorial of the number inside the 
brackets</DIV>
<DIV> </DIV>
<DIV>topo2=topo2.*filter;   %the filter is applied to the second term 
of the series</DIV>
<DIV> </DIV>
<DIV>topo2=constant-topo2; %the new topography approach is computed substracting 
the first term to the second one</DIV>
<DIV> </DIV>
<DIV>topoinverse2=real(ifft2(topo2));   % the real part of the 2-D 
inverse FFT provides the new topography approach in space domain</DIV>
<DIV> </DIV>
<DIV>diference2=topoinverse2-topoinverse;  % this compute the diference 
between the two conscutive topography 
approaches<BR>diference2=diference2.^2;<BR>rms2=sqrt(sum(sum(diference2))/(2*(numrows*numcolumns)));   
%it computes the rms error between the last two iterations<BR>iter=2;  %the 
number of the iteration reached is stored in variable 
'iter'<BR>rms=rms2;<BR>finaltopoinverse=topoinverse2;  %the new topography 
approach is stored in the variable 'finaltpoinverse'<BR> if 
rms2>=criterio   %it determines if the rms error is greater than 
the convergence criterion<BR>    
topo3=(((frequencytotal.^(1))./(prod(1:2))).*(fft2(finaltopoinverse.^2)))+((((frequencytotal.^(2))./(prod(1:3))).*(fft2(finaltopoinverse.^3))));   
%the third term of the series is computed using the new topography 
approach</DIV>
<DIV> </DIV>
<DIV><BR>topo3=topo3.*filter; %the filter is applied to the third term of the 
series</DIV>
<DIV> </DIV>
<DIV> topo3=constant-topo3;  %the new topography approach is 
computed<BR>      
topoinverse3=real(ifft2(topo3));  %and transformed to space 
domain<BR>      
diference3=topoinverse3-topoinverse2;  %the diference between the last two 
topography approaches is computed<BR>      
diference3=diference3.^2;<BR>      
rms3=sqrt(sum(sum(diference3))/(2*(numrows*numcolumns)));  %the new rms 
error is computed<BR>      
rms=rms3;<BR>      finaltopoinverse=topoinverse3;  
%the new topography approach is stored in variable 
finaltopoinverse<BR>      iter=3;  %it indicates 
the iteration step reached<BR>  if rms3>=criterio  %it determines 
if the convergence criterion has been reached or not<BR>    
 topo4=(((frequencytotal.^(1))/(prod(1:2))).*(fft2(topoinverse3.^2)))+((((frequencytotal.^(2))/(prod(1:3))).*(fft2(topoinverse3.^3))))+((((frequencytotal.^(3))/(prod(1:4))).*(fft2(topoinverse3.^4))));  
%the fourth term of the series is computed, and so on...</DIV>
<DIV> </DIV>
<DIV>topo4=topo4.*filter;</DIV>
<DIV> </DIV>
<DIV>     topo4=constant-topo4;<BR>    
 topoinverse4=real(ifft2(topo4));<BR>    
 diference4=topoinverse4-topoinverse3;<BR>    
 diference4=diference4.^2;<BR>      
rms4=sqrt(sum(sum(diference4))/(2*(numrows*numcolumns)));<BR>      
rms=rms4;<BR>      
finaltopoinverse=topoinverse4;<BR>      
iter=4;<BR>       if 
rms4>=criterio<BR>         
topo5=(((frequencytotal.^(1))/(prod(1:2))).*(fft2(topoinverse4.^2)))+((((frequencytotal.^(2))/(prod(1:3))).*(fft2(topoinverse4.^3))))+((((frequencytotal.^(3))/(prod(1:4))).*(fft2(topoinverse4.^4))))+((((frequencytotal.^(4))/(prod(1:5))).*(fft2(topoinverse4.^5))));</DIV>
<DIV> </DIV>
<DIV>topo5=topo5.*filter;</DIV>
<DIV> </DIV>
<DIV>      topo5=constant-topo5;<BR>    
  topoinverse5=real(ifft2(topo5));<BR>    
  diference5=topoinverse5-topoinverse4;<BR>    
  diference5=diference5.^2;<BR>      
 rms5=sqrt(sum(sum(diference5))/(2*(numrows*numcolumns)));<BR>         
rms=rms5;<BR>         
finaltopoinverse=topoinverse5;<BR>         
iter=5;<BR>         if 
rms5>=criterio<BR>            
topo6=(((frequencytotal.^(1))/(prod(1:2))).*(fft2(topoinverse5.^2)))+((((frequencytotal.^(2))/(prod(1:3))).*(fft2(topoinverse5.^3))))+((((frequencytotal.^(3))/(prod(1:4))).*(fft2(topoinverse5.^4))))+((((frequencytotal.^(4))/(prod(1:5))).*(fft2(topoinverse5.^5))))+((((frequencytotal.^(5))/(prod(1:6))).*(fft2(topoinverse5.^6))));</DIV>
<DIV> </DIV>
<DIV>topo6=topo6.*filter;</DIV>
<DIV> </DIV>
<DIV>    
   topo6=constant-topo6;<BR>    
   topoinverse6=real(ifft2(topo6));<BR>    
   diference6=topoinverse6-topoinverse5;<BR>    
   diference6=diference6.^2;<BR>      
  rms6=sqrt(sum(sum(diference6))/(2*(numrows*numcolumns)));<BR>            
rms=rms6;<BR>            
finaltopoinverse=topoinverse6;<BR>            
iter=6;<BR>            if 
rms6>=criterio<BR>               
topo7=(((frequencytotal.^(1))/(prod(1:2))).*(fft2(topoinverse6.^2)))+((((frequencytotal.^(2))/(prod(1:3))).*(fft2(topoinverse6.^3))))+((((frequencytotal.^(3))/(prod(1:4))).*(fft2(topoinverse6.^4))))+((((frequencytotal.^(4))/(prod(1:5))).*(fft2(topoinverse6.^5))))+((((frequencytotal.^(5))/(prod(1:6))).*(fft2(topoinverse6.^6))))+((((frequencytotal.^(6))/(prod(1:7))).*(fft2(topoinverse6.^7))));</DIV>
<DIV> </DIV>
<DIV>topo7=topo7.*filter;</DIV>
<DIV> </DIV>
<DIV>    
    topo7=constant-topo7;<BR>    
    topoinverse7=real(ifft2(topo7));<BR>    
    diference7=topoinverse7-topoinverse6;<BR>    
    diference7=diference7.^2;<BR>      
   rms7=sqrt(sum(sum(diference7))/(2*(numrows*numcolumns)));<BR>               
rms=rms7;<BR>               
finaltopoinverse=topoinverse7;<BR>               
iter=7;<BR>               
if 
rms7>=criterio<BR>                  
topo8=(((frequencytotal.^(1))/(prod(1:2))).*(fft2(topoinverse7.^2)))+((((frequencytotal.^(2))/(prod(1:3))).*(fft2(topoinverse7.^3))))+((((frequencytotal.^(3))/(prod(1:4))).*(fft2(topoinverse7.^4))))+((((frequencytotal.^(4))/(prod(1:5))).*(fft2(topoinverse7.^5))))+((((frequencytotal.^(5))/(prod(1:6))).*(fft2(topoinverse7.^6))))+((((frequencytotal.^(6))/(prod(1:7))).*(fft2(topoinverse7.^7))))+((((frequencytotal.^(7))/(prod(1:8))).*(fft2(topoinverse7.^8))));</DIV>
<DIV> </DIV>
<DIV>topo8=topo8.*filter;<BR>    
     topo8=constant-topo8;<BR>    
     topoinverse8=real(ifft2(topo8));<BR>    
     diference8=topoinverse8-topoinverse7;<BR>    
     diference8=diference8.^2;<BR>      
    rms8=sqrt(sum(sum(diference8))/(2*(numrows*numcolumns)));<BR>                  
rms=rms8;<BR>                  
finaltopoinverse=topoinverse8;<BR>                  
iter=8;<BR>                  
if 
rms8>=criterio<BR>                     
topo9=(((frequencytotal.^(1))/(prod(1:2))).*(fft2(topoinverse8.^2)))+((((frequencytotal.^(2))/(prod(1:3))).*(fft2(topoinverse8.^3))))+((((frequencytotal.^(3))/(prod(1:4))).*(fft2(topoinverse8.^4))))+((((frequencytotal.^(4))/(prod(1:5))).*(fft2(topoinverse8.^5))))+((((frequencytotal.^(5))/(prod(1:6))).*(fft2(topoinverse8.^6))))+((((frequencytotal.^(6))/(prod(1:7))).*(fft2(topoinverse8.^7))))+((((frequencytotal.^(7))/(prod(1:8))).*(fft2(topoinverse8.^8))))+((((frequencytotal.^(8))/(prod(1:9))).*(fft2(topoinverse8.^9))));</DIV>
<DIV> </DIV>
<DIV>topo9=topo9.*filter;<BR>    
      topo9=constant-topo9;<BR>    
      topoinverse9=real(ifft2(topo9));<BR>    
      diference9=topoinverse9-topoinverse8;<BR>    
      diference9=diference9.^2;<BR>      
     rms9=sqrt(sum(sum(diference9))/(2*(numrows*numcolumns)));<BR>                     
rms=rms9;<BR>                     
finaltopoinverse=topoinverse9;<BR>                     
iter=9;<BR>                     
if 
rms9>=criterio<BR>                        
topo10=(((frequencytotal.^(1))/(prod(1:2))).*(fft2(topoinverse9.^2)))+((((frequencytotal.^(2))/(prod(1:3))).*(fft2(topoinverse9.^3))))+((((frequencytotal.^(3))/(prod(1:4))).*(fft2(topoinverse9.^4))))+((((frequencytotal.^(4))/(prod(1:5))).*(fft2(topoinverse9.^5))))+((((frequencytotal.^(5))/(prod(1:6))).*(fft2(topoinverse9.^6))))+((((frequencytotal.^(6))/(prod(1:7))).*(fft2(topoinverse9.^7))))+((((frequencytotal.^(7))/(prod(1:8))).*(fft2(topoinverse9.^8))))+((((frequencytotal.^(8))/(prod(1:9))).*(fft2(topoinverse9.^9))))+((((frequencytotal.^(9))/(prod(1:10))).*(fft2(topoinverse9.^10))));</DIV>
<DIV> </DIV>
<DIV>topo10=topo10.*filter;<BR>    
      topo10=constant-topo10;<BR>    
      topoinverse10=real(ifft2(topo10));<BR>    
      diference10=topoinverse10-topoinverse9;<BR>    
      diference10=diference10.^2;<BR>                     
rms10=sqrt(sum(sum(diference10))/(2*(numrows*numcolumns)));<BR>                     
rms=rms10;<BR>                     
finaltopoinverse=topoinverse10;<BR>                     
iter=10;<BR>                   
end<BR>                
end<BR>             
end<BR>          
end<BR>       end<BR>    
end<BR> end<BR>end;</DIV>
<DIV> </DIV>
<DIV>iter     %it displays the iteration number at which the 
process has stopped<BR>rms      %it displays the rms 
error at the end of the iterative process <BR>figure   %it opens a new 
figure<BR>surf(finaltopoinverse+z0)   %it draws the final topography 
relief map obtained adding the mean reference depth<BR>title('Topography of the 
inverted interface obtained from the Bouguer gravity map (Km)') %this is the 
title of the new graph<BR>%the matrix is transposed to their output to a surfer 
file<BR>trans=(finaltopoinverse+z0)';<BR>fid=fopen(topoout, 
'w');<BR>fprintf(fid,'%6.2f\n',trans);  %this opens and writes the output 
matrix in a ASCII file with a single data column with 
depths<BR>fclose(fid);</DIV>
<DIV> </DIV>
<DIV>%At this point, the forward modelling starts. The final topography obtained 
before is used in the Parker's formula in order to compute the gravity anomaly 
due to that interface. The number of terms of the series is the same that the 
number of iterations reached in the inverse modelling.</DIV>
<DIV> </DIV>
<DIV>sumas=1; %initialize variables. The maximum number of sums is 10, because 
the maximum number of iterations in the inverse modelling is 
10.<BR>sumtotal=[];<BR>sum2=[];<BR>sum3=[];<BR>sum4=[];<BR>sum5=[];<BR>sum6=[];<BR>sum7=[];<BR>sum8=[];<BR>sum9=[];<BR>sum10=[];<BR>constantbouinv=-(2*pi*6.67*contrast)*(exp(-(z0*100000).*(frequencytotal.*(1/100000))));  
%the constant term of the series is computed<BR>%The sum of fourier transforms 
starts here. At each step, a new term of the series is computed and added to the 
previous one. Then, if the number of sums is lower than the number of 
iterations, the process continues.</DIV>
<DIV> </DIV>
<DIV>sumtotal=(((frequencytotal.^(0))/(prod(1:1))).*(fft2(finaltopoinverse.^1)));<BR>sums=2;<BR>if 
sums<=iter<BR>   
sum2=(((frequencytotal.^(1))/(prod(1:2))).*(fft2(finaltopoinverse.^2)));<BR>   
sums=3;<BR>   sumtotal=sumtotal+sum2;<BR>   if 
sumas<=iter<BR>      
sum3=((((frequencytotal.^(2))/(prod(1:3))).*(fft2(finaltopoinverse.^3))));<BR>      
sums=4;<BR>      
sumtotal=sumtotal+sum3;<BR>      if 
sums<=iter<BR>         
sum4=((((frequencytotal.^(3))/(prod(1:4))).*(fft2(finaltopoinverse.^4))));<BR>         
sums=5;<BR>         
sumtotal=sumtotal+sum4;<BR>         if 
sums<=iter<BR>            
sum5=((((frequencytotal.^(4))/(prod(1:5))).*(fft2(finaltopoinverse.^5))));    
  <BR>            
sums=6;<BR>            
sumtotal=sumtotal+sum5;<BR>            
if 
sums<=iter<BR>               
sum6=((((frequencytotal.^(5))/(prod(1:6))).*(fft2(finaltopoinverse.^6))));<BR>               
sums=7;<BR>               
sumtotal=sumtotal+sum6;<BR>               
if 
sums<=iter<BR>                  
sum7=((((frequencytotal.^(6))/(prod(1:7))).*(fft2(finaltopoinverse.^7))));<BR>                  
sums=8;<BR>                  
sumtotal=sumtotal+sum7;<BR>                  
if 
sums<=iter<BR>                     
sum8=((((frequencytotal.^(7))/(prod(1:8))).*(fft2(finaltopoinverse.^8))));<BR>                     
sums=9;<BR>                     
sumtotal=sumtotal+sum8;<BR>                     
if 
sums<=iter<BR>                        
sum9=((((frequencytotal.^(8))/(prod(1:9))).*(fft2(finaltopoinverse.^9))));<BR>                        
sums=10;<BR>                        
sumtotal=sumtotal+sum9;<BR>                        
if 
sums<=iter<BR>                           
sum10=((((frequencytotal.^(9))/(prod(1:10))).*(fft2(finaltopoinverse.^10))));<BR>                           
sumtotal=sumtotal+sum10;<BR>                        
end<BR>                     
end<BR>                  
end<BR>               
end<BR>            
end<BR>         
end<BR>      end<BR>   end<BR>end</DIV>
<DIV> </DIV>
<DIV><BR>sumtotal=sumtotal.*constantbouinv;  %after the summatory, the 
final map in frequency domain is computed multipliying the constant 
term<BR>sumtotal(1,1)=0;<BR>bouinv=real(ifft2(sumtotal));  %the inverse 2-D 
FFT provides the gravity map in space domain</DIV>
<DIV> </DIV>
<DIV>figure %opens a new figure<BR>surf(bouinv)  %it displays the gravity 
map obtained only for visualization<BR>title('Gravity map due to the interface 
calculated (mGal)') %this is the title of the new figure<BR>figure  %opens 
a new figure<BR>surf(bouinv-mapabou) %it displays the map with the diference 
between the gravity map input and the gravity map of the forward 
modelling<BR>title('Diference between the input gravity map and the one due to 
the calculated interface (mGal)')  %this is the title of the new 
figure</DIV>
<DIV> </DIV>
<DIV>bouinvtras=bouinv';  %the matrix is transposed to their output to a 
surfer file<BR>fid=fopen(bouinverted, 
'w');<BR>fprintf(fid,'%6.2f\n',bouinvtras+meangravity);  %this opens and 
writes the output matrix (adding the mean gravity value substracted prior to the 
FFT) in a ASCII file with a single data column with gravity 
anomaly<BR>fclose(fid);</DIV></FONT></BODY></HTML>