[Scilab-users] Matlab to Scilab conversion issues
arctica1963
arctica1963 at gmail.com
Sun Sep 16 23:09:09 CEST 2012
Hello,
I have tried to run the mfile2sci conversion on the Matlab file with a few
compatibility issues that I isolated from the converted file.
tukeywin - a tapered cosine window which was a function in a Matlab toolbox
[ cannot see an obvious Scilab function(s) ]
ifft2 - 2D inverse discrete fourier transform [ is this the same as Scilab
ifft ? ]
figure - creates a graphics object window [ a Scilab equivalent to open a
new graphics window? ]
The original filename began with a number so had to change that to get the
conversion to run. I have added the converted file below:
------------------------------------------------------------------------------------------------
// Display mode
mode(0);
// Display warning for floating point exception
ieee(1);
//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 = mtlb_fopen(bouin);
// L.40: No simple equivalent, so mtlb_fscanf() is called.
bou = mtlb_fscanf(fid,"%f",[numcolumns,numrows]);
//open the input gravity file, read it and stores it in variable ''bou''
mclose(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
// !! L.52: Matlab toolbox(es) function tukeywin not converted, original
calling sequence used
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''
// !! L.53: Matlab toolbox(es) function tukeywin not converted, original
calling sequence used
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 =
mtlb_double(mtlb_e(wrows,":"))*mtlb_double(mtlb_0(mtlb_e(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
// !! L.64: Matlab function figure not yet converted.
// !! L.64: Matlab function figure not yet converted, original calling
sequence used.
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 = frequency(:,$:-1:1);
frequency3 = frequency($:-1:1,:);
%v0 = frequency($:-1:1,:);
frequency4 = %v0(:,$:-1:1);
entero = round(numcolumns/2);
if (numcolumns/2-entero)==0 then
frequency2(:,1) = [];
frequency3(1,:) = [];
frequency4(:,1) = [];
frequency4(1,:) = [];
frequencytotal = [frequency,frequency2;frequency3,frequency4];
else
frequencytotal = [frequency,frequency2;frequency3,frequency4];
end;
frequencytotal($,:) = [];
frequencytotal(:,$) = [];
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 then
filter(f,g) = 1;
elseif frequencytotal(f,g)<SH then
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
// !! L.114: Matlab function ifft2 not yet converted, original calling
sequence used.
topoinverse = real(mtlb_double(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 = mtlb_s(constant,topo2);
//the new topography approach is computed substracting the first term to the
second one
// !! L.125: Matlab function ifft2 not yet converted, original calling
sequence used.
topoinverse2 = real(mtlb_double(ifft2(topo2)));
// the real part of the 2-D inverse FFT provides the new topography approach
in space domain
diference2 = mtlb_s(topoinverse2,topoinverse);
// this compute the diference between the two conscutive topography
approaches
diference2 = diference2 .^2;
rms2 = sqrt(mtlb_sum(mtlb_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 mtlb_logic(rms2,">=",criterio) then
//it determines if the rms error is greater than the convergence criterion
topo3 = mtlb_a(((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 = mtlb_s(constant,topo3);
//the new topography approach is computed
// !! L.140: Matlab function ifft2 not yet converted, original calling
sequence used.
topoinverse3 = real(mtlb_double(ifft2(topo3)));
//and transformed to space domain
diference3 = mtlb_s(topoinverse3,topoinverse2);
//the diference between the last two topography approaches is computed
diference3 = diference3 .^2;
rms3 = sqrt(mtlb_sum(mtlb_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 mtlb_logic(rms3,">=",criterio) then //it determines if the convergence
criterion has been reached or not
topo4 = mtlb_a(mtlb_a(((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 = mtlb_s(constant,topo4);
// !! L.153: Matlab function ifft2 not yet converted, original calling
sequence used.
topoinverse4 = real(mtlb_double(ifft2(topo4)));
diference4 = mtlb_s(topoinverse4,topoinverse3);
diference4 = diference4 .^2;
rms4 = sqrt(mtlb_sum(mtlb_sum(diference4))/(2*(numrows*numcolumns)));
rms = rms4;
finaltopoinverse = topoinverse4;
iter = 4;
if mtlb_logic(rms4,">=",criterio) then
topo5 = mtlb_a(mtlb_a(mtlb_a(((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 = mtlb_s(constant,topo5);
// !! L.166: Matlab function ifft2 not yet converted, original calling
sequence used.
topoinverse5 = real(mtlb_double(ifft2(topo5)));
diference5 = mtlb_s(topoinverse5,topoinverse4);
diference5 = diference5 .^2;
rms5 = sqrt(mtlb_sum(mtlb_sum(diference5))/(2*(numrows*numcolumns)));
rms = rms5;
finaltopoinverse = topoinverse5;
iter = 5;
if mtlb_logic(rms5,">=",criterio) then
topo6 = mtlb_a(mtlb_a(mtlb_a(mtlb_a(((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 = mtlb_s(constant,topo6);
// !! L.179: Matlab function ifft2 not yet converted, original
calling sequence used.
topoinverse6 = real(mtlb_double(ifft2(topo6)));
diference6 = mtlb_s(topoinverse6,topoinverse5);
diference6 = diference6 .^2;
rms6 = sqrt(mtlb_sum(mtlb_sum(diference6))/(2*(numrows*numcolumns)));
rms = rms6;
finaltopoinverse = topoinverse6;
iter = 6;
if mtlb_logic(rms6,">=",criterio) then
topo7 = mtlb_a(mtlb_a(mtlb_a(mtlb_a(mtlb_a(((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 = mtlb_s(constant,topo7);
// !! L.192: Matlab function ifft2 not yet converted, original
calling sequence used.
topoinverse7 = real(mtlb_double(ifft2(topo7)));
diference7 = mtlb_s(topoinverse7,topoinverse6);
diference7 = diference7 .^2;
rms7 =
sqrt(mtlb_sum(mtlb_sum(diference7))/(2*(numrows*numcolumns)));
rms = rms7;
finaltopoinverse = topoinverse7;
iter = 7;
if mtlb_logic(rms7,">=",criterio) then
topo8 =
mtlb_a(mtlb_a(mtlb_a(mtlb_a(mtlb_a(mtlb_a(((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 = mtlb_s(constant,topo8);
// !! L.204: Matlab function ifft2 not yet converted, original
calling sequence used.
topoinverse8 = real(mtlb_double(ifft2(topo8)));
diference8 = mtlb_s(topoinverse8,topoinverse7);
diference8 = diference8 .^2;
rms8 =
sqrt(mtlb_sum(mtlb_sum(diference8))/(2*(numrows*numcolumns)));
rms = rms8;
finaltopoinverse = topoinverse8;
iter = 8;
if mtlb_logic(rms8,">=",criterio) then
topo9 =
mtlb_a(mtlb_a(mtlb_a(mtlb_a(mtlb_a(mtlb_a(mtlb_a(((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 = mtlb_s(constant,topo9);
// !! L.216: Matlab function ifft2 not yet converted, original
calling sequence used.
topoinverse9 = real(mtlb_double(ifft2(topo9)));
diference9 = mtlb_s(topoinverse9,topoinverse8);
diference9 = diference9 .^2;
rms9 =
sqrt(mtlb_sum(mtlb_sum(diference9))/(2*(numrows*numcolumns)));
rms = rms9;
finaltopoinverse = topoinverse9;
iter = 9;
if mtlb_logic(rms9,">=",criterio) then
topo10 =
mtlb_a(mtlb_a(mtlb_a(mtlb_a(mtlb_a(mtlb_a(mtlb_a(mtlb_a(((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 = mtlb_s(constant,topo10);
// !! L.228: Matlab function ifft2 not yet converted,
original calling sequence used.
topoinverse10 = real(mtlb_double(ifft2(topo10)));
diference10 = mtlb_s(topoinverse10,topoinverse9);
diference10 = diference10 .^2;
rms10 =
sqrt(mtlb_sum(mtlb_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
// !! L.246: Matlab function figure not yet converted.
// !! L.246: Matlab function figure not yet converted, original calling
sequence used.
figure;
//it opens a new figure
surf(mtlb_a(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 = mtlb_a(finaltopoinverse,z0)';
fid = mtlb_fopen(topoout,"w");
// L.252: No simple equivalent, so mtlb_fprintf() is called.
mtlb_fprintf(fid,"%6.2f\n",trans);
//this opens and writes the output matrix in a ASCII file with a single data
column with depths
mclose(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,firstnonsingleton(1:1)))
.*fft2(finaltopoinverse .^1);
sums = 2;
if sums<=iter then
sum2 = ((frequencytotal .^1)/prod(1:2)) .*fft2(finaltopoinverse .^2);
sums = 3;
sumtotal = mtlb_a(sumtotal,sum2);
if sumas<=iter then
sum3 = ((frequencytotal .^2)/prod(1:3)) .*fft2(finaltopoinverse .^3);
sums = 4;
sumtotal = mtlb_a(sumtotal,sum3);
if sums<=iter then
sum4 = ((frequencytotal .^3)/prod(1:4)) .*fft2(finaltopoinverse .^4);
sums = 5;
sumtotal = mtlb_a(sumtotal,sum4);
if sums<=iter then
sum5 = ((frequencytotal .^4)/prod(1:5)) .*fft2(finaltopoinverse
.^5);
sums = 6;
sumtotal = mtlb_a(sumtotal,sum5);
if sums<=iter then
sum6 = ((frequencytotal .^5)/prod(1:6)) .*fft2(finaltopoinverse
.^6);
sums = 7;
sumtotal = mtlb_a(sumtotal,sum6);
if sums<=iter then
sum7 = ((frequencytotal .^6)/prod(1:7)) .*fft2(finaltopoinverse
.^7);
sums = 8;
sumtotal = mtlb_a(sumtotal,sum7);
if sums<=iter then
sum8 = ((frequencytotal .^7)/prod(1:8))
.*fft2(finaltopoinverse .^8);
sums = 9;
sumtotal = mtlb_a(sumtotal,sum8);
if sums<=iter then
sum9 = ((frequencytotal .^8)/prod(1:9))
.*fft2(finaltopoinverse .^9);
sums = 10;
sumtotal = mtlb_a(sumtotal,sum9);
if sums<=iter then
sum10 = ((frequencytotal .^9)/prod(1:10))
.*fft2(finaltopoinverse .^10);
sumtotal = mtlb_a(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;
// !! L.321: Matlab function ifft2 not yet converted, original calling
sequence used.
bouinv = real(mtlb_double(ifft2(sumtotal)));
//the inverse 2-D FFT provides the gravity map in space domain
// !! L.323: Matlab function figure not yet converted.
// !! L.323: Matlab function figure not yet converted, original calling
sequence used.
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
// !! L.326: Matlab function figure not yet converted.
// !! L.326: Matlab function figure not yet converted, original calling
sequence used.
figure;
//opens a new figure
surf(mtlb_s(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 = mtlb_fopen(bouinverted,"w");
// L.332: No simple equivalent, so mtlb_fprintf() is called.
mtlb_fprintf(fid,"%6.2f\n",mtlb_a(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
mclose(fid);
--------------------------------------------------------------------------------
If anyone can provide pointers on the equivalent conversion for Scilab or
other guidance, that would be much appreciated.
I have not used Matlab myself, so not sure of the way to proceed.
Cheers
Lester
--
View this message in context: http://mailinglists.scilab.org/Matlab-to-Scilab-conversion-issues-tp4024825.html
Sent from the Scilab users - Mailing Lists Archives mailing list archive at Nabble.com.
More information about the users
mailing list