Hello,
You just have to replace "x" by "wavelength" and "y" by "ln_power".
The slopes are the first two component of the optimal vector "popt" :
clearclf()// Read data - wavelength (in km)), power, 1 standard
deviation// Unknown data length; 3 columns -default space delimited //
PSD_wavelength.dat from GMT grdfft radially averaged power spectra
data = read('PSD_wavelength.dat',-1,3); wavelength = data(:,1); power
= data(:,2); std_dev1 = data(:,3); ln_power = log(power); wavenumber
= 1./wavelength; f=gcf(); //plot(wavenumber, ln_power)
scatter(wavenumber, ln_power,'marker','.') a=gce().children;
a.mark_mode = "on"a.mark_style = 0a.mark_size_unit =
"point"a.mark_size=3 xlabel ('wavenumber k (km-1)')ylabel ('Log
(Power)') function y=fun(x, param) a = param(1); b = param(2);
theta = param(3); phi = param(4); c =
(x<theta)*a+(x>=theta)*b; y = c.*(x-theta)+phi; endfunction
function r=residual(param, x, y, f) r = sum((y-f(x,param)).^2);
endfunction function [r, dr, ind]=costf(p, ind, x, y, f) r =
residual(p,x,y,f); dr = numderivative(list(residual,x,y,fun),p);
endfunction [ropt,popt] =
optim(list(costf,wavenumber,ln_power,fun),[0,0,mean(wavenumber),mean(ln_power)]);
plot(wavenumber,fun(wavenumber,popt),'r',popt(3),popt(4),'xr')
arctica1963 <arctica1...@gmail.com> a écrit :
Hello,
Thanks for the idea and suggestions. Not too sure how to apply it, if you
could give some pointers on the attached data and code. The ultimate idea is
to get the slopes of the straight line segments. Many thanks, Lester
clear
clf()
// Read data - wavelength (in km)), power, 1 standard deviation
// Unknown data length; 3 columns -default space delimited
// PSD_wavelength.dat from GMT grdfft radially averaged power spectra
data = read('PSD_wavelength.dat',-1,3);
wavelength = data(:,1);
power = data(:,2);
std_dev1 = data(:,3);
ln_power = log(power);
wavenumber = 1./wavelength;
f=gcf();
//plot(wavenumber, ln_power)
scatter(wavenumber, ln_power,'marker','.')
a=gce().children;
a.mark_mode = "on"
a.mark_style = 0
a.mark_size_unit = "point"
a.mark_size=3
xlabel ('wavenumber k (km-1)')
ylabel ('Log (Power)')
PSD_wavelength.dat
<https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/mailinglists.scilab.org/file/t495709/PSD_wavelength.dat>
--
Sent from:
https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
_______________________________________________
users mailing list
users@lists.scilab.orghttps://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/lists.scilab.org/mailman/listinfo/users
_______________________________________________
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users
