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='caf2_maggio10_1725';
str=['load ' file '.dat'];eval(str);
str=['SHGx=' file '(:,1);'];eval(str);
str=['SHGy=' file '(:,2);'];eval(str);                
t=(-10000:0.1:10000);% fs



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

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

%w=(2.122:5e-5*2*pi:2.75);% fs^-1
w=(2.122:5e-5*2*pi:2.75);% 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=3300;   %micron silica e CaF2

% %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;
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;
    i
end
toc

figure(2)
plot(t,real(y1))
% i termini y1^4 sono solo additivi, non si considerano, questo genera
% l'"artificio" di una funzione risposta negativa.
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+1:dimt*3/2+1);
% filtro
A=fft(SHG1);
cut=4000;
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',t,exp(-4*log(2)*t.^2/(2^0.5*Dt)^2),'b',SHGx,SHGy,'r');xlim([-400 400])
% plot(t,SHG1/max(SHG1),'r',t,SHG1filter/max(SHG1filter),'k');
if 0, 
    AA=[t' (yy-Ybkg)']; str=['save ' filecommon '.dat AA -ascii;'];eval(str);
    AA=[t' (SHG1/max(SHG1))' (SHG1filter/max(SHG1filter))']; str=['save shg_simul.dat AA -ascii;'];eval(str);
end;