clear all
file='spettrolaserG3';
str=['load ' file '.txt'];eval(str);
str=['Ls=1e-3*' file '(:,1);'];eval(str);
str=['spec=' file '(:,2);'];eval(str);

file='ca2f_maggio10_1725';
str=['load ' file '.dat'];eval(str);
str=['SHGx=' file '(:,1);'];eval(str);
str=['SHGy=' file '(:,2);'];eval(str);                
t=(-1000:0.1:1000);% fs



Dt=13; % fs
Dw=0.441/Dt*2*pi;

c=0.3;    %micron/fs
lambda0=0.795; %lunghezza d'onda centrale in micron
w0=2*pi*c/lambda0;%pulsazione centrale in fs^-1
kk=1;

w=(w0-w0/kk:1e-4*2*pi:w0+w0/kk);% fs^-1

nn=find(Ls<0.6);spec=spec-mean(spec(nn));spec=spec/max(spec);
specI=interp1(c./Ls*2*pi,spec,w,'linear',0);
specI2=exp(-4*log(2)*(w-w0).^2/Dw^2);

figure(1)
plot(w,specI,'.',w,specI2,'.')


L=0;   %micron

% %silica
n0=1.453;%indice di rifrazione
d1nl=-0.01728;%derivata prima di n in lambda micron^-1
d2nl=0.03988; %derivata seconda di n in lambda micron^-2
d3nl=-0.2388; %derivata terza di n in lambda micron^-3

% % BK7
% n0=1.511;%indice di rifrazione
% d1nl=-0.01456;%derivata prima di n in lambda micron^-1
% d2nl=0.05589; %derivata seconda di n in lambda micron^-2
% d3nl=-0.2877; %derivata terza di n in lambda micron^-3

% % CaF2
% n0=1.4306;
% d1nl=-0.0104;
% d2nl=0.0304;
% d3nl=-0.1662;

% % Calcite ordinario
% n0=1.6488;
% d1nl=-0.031048300320746;
% d2nl=0.081657209660690;
% d3nl=-0.476683581496218;

% % Calcite straordinario
% n0=1.481916136599520;
% d1nl=-0.012620442731315;
% d2nl=0.044786681918652;
% d3nl=-0.232969865621158;
% 
% % Benzene
% n0=1.485267310086479;
% d1nl=-0.029582350419277;
% d2nl=0.115024263003249;
% d3nl=-0.609046368893473;




d1nw=-lambda0^2/2/pi/c*d1nl; %derivata prima indice di rifrazione in funzione di w fs^-1
d2nw=lambda0^3/4/pi^2/c^2*(2*d1nl+lambda0*d2nl); %derivata seconta indice di rifrazione in funzione di w fs^-2
d3nw=-lambda0^4/8/pi^3/c^3*(6*d1nl+6*lambda0*d2nl+lambda0^2*d3nl); %derivata terza indice di rifrazione in funzione di w fs^-3


phi0=n0*w0*L/c;
phi1=L/c*(n0+w0*d1nw);
phi2=L/c*(2*d1nw+w0*d2nw);
phi3=L/c*(3*d2nw+w0*d3nw);

dim=length(w);

A=sqrt(specI);

y1=t*0;
AA=w*0;
phi=w*0;

tic
for i=1:1:dim 
    phi(i)=0*phi0+0*phi1*(w(i)-w0)+phi2*(w(i)-w0).^2/2+phi3*(w(i)-w0).^3/6;
    y1=A(i)*cos(w(i)*t-phi(i))+y1;
    AA(i)=A(i);
    i
end
toc
if 0
    pp=find(t>-300 & t<300);ss=y1(pp);
    nn=find(diff(ss)<=0);ss(nn)=0;
    mm=find(diff(ss)>0);ss(mm)=1;
    gg=find(diff(ss)>0);jj=diff(t(gg));
    plot(jj*2)
    
end
figure(2)
plot(t,y1)
shg1=4*xcorr(y1.^2,y1.^2)+1*(4*xcorr(y1.^3,y1)+4*xcorr(y1,y1.^3));
dimt=length(t)-1;
SHG1=shg1(dimt/2:dimt*3/2);
% figure(2)
% plot(t,SHG1/max(SHG1),t,exp(-4*log(2)*t.^2/(2^0.5*Dt)^2))
% filtro
A=fft(SHG1);
cut=400;
A(cut:end-cut)=SHG1(cut:end-cut)*0;
SHG1filter=real(ifft(A));
figure(3)
% plot(t,SHG1/max(SHG1),t,SHG1filter/max(SHG1filter),SHGx,SHGy)
% plot(t,SHG1/max(SHG1)*4.5,t,SHG1filter/max(SHG1filter),'k',SHGx-0.1,SHGy,'r');xlim([-400 400])
plot(t,SHG1filter/max(SHG1filter),'k',SHGx,SHGy,'r');xlim([-400 400])