#! esempi di manipolazione degli spettri
## necessita del pacchetto signal
## pkg install -forge general control signal
## usiamo anche findpeaksG.m da
## http://terpconnect.umd.edu/~toh/spectrum/PeakFindingandMeasurement.htm
1;

xxx=csvread("CO.csv",2,0);  # leggi lo spettro di CO, in formato CSV, salta 2 righe e 0 colonne

x=flipdim(xxx(:,1));
y=flipdim(xxx(:,2));

ixa=lookup(x,1900);
ixb=lookup(x,2300);

## restringi l'intervallo d'interesse alla banda fondamentale

x1=x(ixa:ixb);
y1=y(ixa:ixb);

## trova i picchi, il segnale dovrebbe essere positivo
y1=y1-min(y);

[~,ixp]=findpeaks(y1,"minpeakheight",0.1);

## meglio, per eliminare i satellitini
[~,ixp]=findpeaks(y1,"minpeakheight",0.1,"minpeakdistance",20);

plot(x1,y1,x1(ixp),y1(ixp),"or","markersize",4)

ixc=lookup(x,4100);
ixd=lookup(x,4400);

x2=x(ixc:ixd);
y2=y(ixc:ixd)-min(y);

[~,ixp2]=findpeaks(y2,"minpeakheight",0.02,"minpeakdistance",20);

plot(x2,y2,x2(ixp2),y2(ixp2),"or","markersize",4)


xxx=csvread("cella.csv",2,0);

xc=flipdim(xxx(:,1));
yc=flipdim(xxx(:,2));

ixca=lookup(xc,3000);
ixcb=lookup(xc,4000);

xc=xc(ixca:ixcb);
yc=yc(ixca:ixcb)-min(yc);

[~,ixcp]=findpeaks(yc,"minpeakheight",0.02);
pc=findpeaksG(xc,yc,0,0.01,100,100)

plot(xc,yc,pc(:,2),pc(:,3),"or")


##------peak fit position

ixp=sort(ixp);
px=x1(ixp);
py=y1(ixp);
pxc=lookup(px,2141);

plot(x1,y1,px(1:pxc),py(1:pxc),"or",px(pxc+1:end),py(pxc+1:end),"ob")

fit_x=[-pxc:-1,1:numel(px)-pxc]';
fit_y=[px(1:pxc);px(pxc+1:end)];

plot(fit_x,fit_y,"xr","markersize",2)
grid

[fit_p1,fit_s1]=polyfit(fit_x,fit_y,1);

plot(fit_x,fit_y,"xr","markersize",2,fit_x,fit_s1.yf)

[fit_p3,fit_s3]=polyfit(fit_x,fit_y,3);

plot(fit_x,fit_y,"xr","markersize",2,fit_x,fit_s1.yf,fit_x,fit_s3.yf)

##---------peak fit shape

rosso=csvread("Rosso.txt",2,0)';

rosso=rosso(:,30:180); # taglio a 600 nm
light=2.99792458e8 #m/s
rosso(1,:)=light/1000 ./ rosso(1,:); #THz

rx=rosso(1,:);
ry=rosso(2,:);

inparms=[0.2;570;20;0.1;610;20];

plot_results(rx,ry,inparms);

fitgaussians(rx,ry,inparms)

inparms=[0.2;570;20;0.1;610;20;0.05;650;20];

plot_results(rx,ry,inparms);

fitgaussians(rx,ry,inparms)

inparms=[0.2;570;25;0.1;610;25;0.05;650;25; 0.01; 690;25];

plot_results(rx,ry,inparms);

fitgaussians(rx,ry,inparms)

inparms=[0.2;572;34;  0.13;612;34;  0.04;650;34;  0.012; 689;34; ];

plot_results(rx,ry,inparms);

parms=fitgaussians(rx,ry,inparms)

inparms=[parms;  0.005;530;36];

plot_results(rx,ry,inparms);

parms=fitgaussians(rx,ry,inparms)

