Convert matlab m file to scilab sci format
Lester Anderson
lester at arctica1.wanadoo.co.uk
Thu Jun 26 19:22:32 CEST 2008
Hello
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
I am not a Matlab user so cannot say why it is failing:
I get the following error:
Startup execution:
loading initial environment
-->mfile2sci();
****** Beginning of mfile2sci() session ******
File to convert: D:/3D_Inver/3dinver.m
Result file path: D:/Temp/
Recursive mode: OFF
Only double values used in M-file: NO
Verbose mode: 3
Generate formated code: NO
M-file reading...
M-file reading: Done
Syntax modification...
Syntax modification: Done
!--error 37
incorrect function at line 1
at line 190 of function mfile2sci called by :
line 64 of function m2sci_gui called by :
line 16 of function mfile2sci called by :
mfile2sci();
-->
I have listed the file here in case someone can see what the issue is
Thanks
Lester
%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=112; %e. g. 112 data in Y direction
numcolumns=112; % e. g. 112 data in X direction
longx=666; %e.g. 666 Km data length in X direction
longy=666; %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.01; %smaller cut-off frequency (1/wavelength in Km) e. g. 0.01, i.e. 1/100 Km
SH=0.012; %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(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
% 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(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'
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'
w2 =wrows(:)* wcolumns(:).'; %this generates a 2-D cosine Tukey window multipliying the 1-D windows
bou=bou.*w2'; %the original gravity input matrix, previously demeaned, is multiplied by the cosine window
mapabou=bou'; %the original gravity input matrix after demeaning is transposed
surf(mapabou) %then it is drawn only for visualization
title('Observed gravity anomaly map (mGal)') %this is the title of the graph
fftbou=fft2(mapabou); %the 2-D FFT of the gravity input matrix is computed after demeaning
spectrum=abs(fftbou(1:numrows/2, 1:numcolumns/2)); %this computes the amplitude spectrum
figure
surf(spectrum) %and the spectrum is drawn only for visualization
title('Amplitude spectrum of the Observed gravity map') %this is the title of the new graph
%the matrix with the frequencies of every harmonic is computed
for f=1:abs((numrows/2)+1);
for g=1:abs((numcolumns/2)+1);
frequency(f, g)=sqrt(((f-1)/longx)^2+((g-1)/longy)^2);
end;
end;
%the matrix of the negative frequencies is also computed
frequency2=fliplr(frequency);
frequency3=flipud(frequency);
frequency4=fliplr(flipud(frequency));
entero=round(numcolumns/2);
if ((numcolumns/2)- entero)==0
frequency2(:,1)=[];
frequency3(1,:)=[];
frequency4(:,1)=[];
frequency4(1,:)=[];
frequencytotal=[frequency frequency2;frequency3 frequency4];
else
frequencytotal=[frequency frequency2;frequency3 frequency4];
end
frequencytotal(end,:)=[];
frequencytotal(:,end)=[];
frequencytotal=frequencytotal.*(2*pi); %the frequency (1/wavelength) matrix is transformed to wavenumber (2*pi/wavelength) matrix
%The iterative process starts here
%The first term of the series, that is constant, is computed and stored in variable 'constant'
up=-(fftbou.*(exp((z0*100000)*(frequencytotal.*(1/100000)))));
down=(2*pi*6.67*contrast);
constant=up./down;
%The high-cut filter is constructed
filter=frequencytotal.*0; %the filter matrix is set to zero
frequencytotal=frequencytotal./(2*pi);
for f=1:numrows;
for g=1:numcolumns;
if frequencytotal(f,g)<WH
filter(f,g)=1;
elseif frequencytotal(f,g)<SH
filter(f,g)=0.5.*(1+cos((((2*pi)*frequencytotal(f,g))-(2*pi*WH))/(2*(SH-WH))));
else
filter(f,g)=0;
end
end;
end;
constant=constant.*filter; %the filter is applied to the first term of the series
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'
frequencytotal=frequencytotal.*(2*pi);
%It starts the computation of the second term of the series with a maximum of 10 iterations
topo2=(((frequencytotal.^(1))./(prod(1:2))).*(fft2(topoinverse.^2))); %the function 'prod' computes the factorial of the number inside the brackets
topo2=topo2.*filter; %the filter is applied to the second term of the series
topo2=constant-topo2; %the new topography approach is computed substracting the first term to the second one
topoinverse2=real(ifft2(topo2)); % the real part of the 2-D inverse FFT provides the new topography approach in space domain
diference2=topoinverse2-topoinverse; % this compute the diference between the two conscutive topography approaches
diference2=diference2.^2;
rms2=sqrt(sum(sum(diference2))/(2*(numrows*numcolumns))); %it computes the rms error between the last two iterations
iter=2; %the number of the iteration reached is stored in variable 'iter'
rms=rms2;
finaltopoinverse=topoinverse2; %the new topography approach is stored in the variable 'finaltpoinverse'
if rms2>=criterio %it determines if the rms error is greater than the convergence criterion
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
topo3=topo3.*filter; %the filter is applied to the third term of the series
topo3=constant-topo3; %the new topography approach is computed
topoinverse3=real(ifft2(topo3)); %and transformed to space domain
diference3=topoinverse3-topoinverse2; %the diference between the last two topography approaches is computed
diference3=diference3.^2;
rms3=sqrt(sum(sum(diference3))/(2*(numrows*numcolumns))); %the new rms error is computed
rms=rms3;
finaltopoinverse=topoinverse3; %the new topography approach is stored in variable finaltopoinverse
iter=3; %it indicates the iteration step reached
if rms3>=criterio %it determines if the convergence criterion has been reached or not
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...
topo4=topo4.*filter;
topo4=constant-topo4;
topoinverse4=real(ifft2(topo4));
diference4=topoinverse4-topoinverse3;
diference4=diference4.^2;
rms4=sqrt(sum(sum(diference4))/(2*(numrows*numcolumns)));
rms=rms4;
finaltopoinverse=topoinverse4;
iter=4;
if rms4>=criterio
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))));
topo5=topo5.*filter;
topo5=constant-topo5;
topoinverse5=real(ifft2(topo5));
diference5=topoinverse5-topoinverse4;
diference5=diference5.^2;
rms5=sqrt(sum(sum(diference5))/(2*(numrows*numcolumns)));
rms=rms5;
finaltopoinverse=topoinverse5;
iter=5;
if rms5>=criterio
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))));
topo6=topo6.*filter;
topo6=constant-topo6;
topoinverse6=real(ifft2(topo6));
diference6=topoinverse6-topoinverse5;
diference6=diference6.^2;
rms6=sqrt(sum(sum(diference6))/(2*(numrows*numcolumns)));
rms=rms6;
finaltopoinverse=topoinverse6;
iter=6;
if rms6>=criterio
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))));
topo7=topo7.*filter;
topo7=constant-topo7;
topoinverse7=real(ifft2(topo7));
diference7=topoinverse7-topoinverse6;
diference7=diference7.^2;
rms7=sqrt(sum(sum(diference7))/(2*(numrows*numcolumns)));
rms=rms7;
finaltopoinverse=topoinverse7;
iter=7;
if rms7>=criterio
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))));
topo8=topo8.*filter;
topo8=constant-topo8;
topoinverse8=real(ifft2(topo8));
diference8=topoinverse8-topoinverse7;
diference8=diference8.^2;
rms8=sqrt(sum(sum(diference8))/(2*(numrows*numcolumns)));
rms=rms8;
finaltopoinverse=topoinverse8;
iter=8;
if rms8>=criterio
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))));
topo9=topo9.*filter;
topo9=constant-topo9;
topoinverse9=real(ifft2(topo9));
diference9=topoinverse9-topoinverse8;
diference9=diference9.^2;
rms9=sqrt(sum(sum(diference9))/(2*(numrows*numcolumns)));
rms=rms9;
finaltopoinverse=topoinverse9;
iter=9;
if rms9>=criterio
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))));
topo10=topo10.*filter;
topo10=constant-topo10;
topoinverse10=real(ifft2(topo10));
diference10=topoinverse10-topoinverse9;
diference10=diference10.^2;
rms10=sqrt(sum(sum(diference10))/(2*(numrows*numcolumns)));
rms=rms10;
finaltopoinverse=topoinverse10;
iter=10;
end
end
end
end
end
end
end
end;
iter %it displays the iteration number at which the process has stopped
rms %it displays the rms error at the end of the iterative process
figure %it opens a new figure
surf(finaltopoinverse+z0) %it draws the final topography relief map obtained adding the mean reference depth
title('Topography of the inverted interface obtained from the Bouguer gravity map (Km)') %this is the title of the new graph
%the matrix is transposed to their output to a surfer file
trans=(finaltopoinverse+z0)';
fid=fopen(topoout, 'w');
fprintf(fid,'%6.2f\n',trans); %this opens and writes the output matrix in a ASCII file with a single data column with depths
fclose(fid);
%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.
sumas=1; %initialize variables. The maximum number of sums is 10, because the maximum number of iterations in the inverse modelling is 10.
sumtotal=[];
sum2=[];
sum3=[];
sum4=[];
sum5=[];
sum6=[];
sum7=[];
sum8=[];
sum9=[];
sum10=[];
constantbouinv=-(2*pi*6.67*contrast)*(exp(-(z0*100000).*(frequencytotal.*(1/100000)))); %the constant term of the series is computed
%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.
sumtotal=(((frequencytotal.^(0))/(prod(1:1))).*(fft2(finaltopoinverse.^1)));
sums=2;
if sums<=iter
sum2=(((frequencytotal.^(1))/(prod(1:2))).*(fft2(finaltopoinverse.^2)));
sums=3;
sumtotal=sumtotal+sum2;
if sumas<=iter
sum3=((((frequencytotal.^(2))/(prod(1:3))).*(fft2(finaltopoinverse.^3))));
sums=4;
sumtotal=sumtotal+sum3;
if sums<=iter
sum4=((((frequencytotal.^(3))/(prod(1:4))).*(fft2(finaltopoinverse.^4))));
sums=5;
sumtotal=sumtotal+sum4;
if sums<=iter
sum5=((((frequencytotal.^(4))/(prod(1:5))).*(fft2(finaltopoinverse.^5))));
sums=6;
sumtotal=sumtotal+sum5;
if sums<=iter
sum6=((((frequencytotal.^(5))/(prod(1:6))).*(fft2(finaltopoinverse.^6))));
sums=7;
sumtotal=sumtotal+sum6;
if sums<=iter
sum7=((((frequencytotal.^(6))/(prod(1:7))).*(fft2(finaltopoinverse.^7))));
sums=8;
sumtotal=sumtotal+sum7;
if sums<=iter
sum8=((((frequencytotal.^(7))/(prod(1:8))).*(fft2(finaltopoinverse.^8))));
sums=9;
sumtotal=sumtotal+sum8;
if sums<=iter
sum9=((((frequencytotal.^(8))/(prod(1:9))).*(fft2(finaltopoinverse.^9))));
sums=10;
sumtotal=sumtotal+sum9;
if sums<=iter
sum10=((((frequencytotal.^(9))/(prod(1:10))).*(fft2(finaltopoinverse.^10))));
sumtotal=sumtotal+sum10;
end
end
end
end
end
end
end
end
end
sumtotal=sumtotal.*constantbouinv; %after the summatory, the final map in frequency domain is computed multipliying the constant term
sumtotal(1,1)=0;
bouinv=real(ifft2(sumtotal)); %the inverse 2-D FFT provides the gravity map in space domain
figure %opens a new figure
surf(bouinv) %it displays the gravity map obtained only for visualization
title('Gravity map due to the interface calculated (mGal)') %this is the title of the new figure
figure %opens a new figure
surf(bouinv-mapabou) %it displays the map with the diference between the gravity map input and the gravity map of the forward modelling
title('Diference between the input gravity map and the one due to the calculated interface (mGal)') %this is the title of the new figure
bouinvtras=bouinv'; %the matrix is transposed to their output to a surfer file
fid=fopen(bouinverted, 'w');
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
fclose(fid);
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20080626/f8ea228e/attachment.htm>
More information about the users
mailing list