clear all
           
t=(-1000:0.1:1000);% fs



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

c=0.3;    %micron/fs
lambda0=0.8; %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


specI2=exp(-4*log(2)*(w-w0).^2/Dw^2);

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


L=500;   %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(specI2);

y1=t*0;
y2=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;
    y2=A(i)*cos(w(i)*t)+y2;
    AA(i)=A(i);
    i
end
toc

figure(2)
plot(t,y1,t,y2)
mac1=xcorr(y1,y1);
mac2=xcorr(y2,y2);

shg1=4*xcorr(y1.^2,y1.^2)+4*xcorr(y1.^3,y1)+4*xcorr(y1,y1.^3);
shg2=4*xcorr(y2.^2,y2.^2)+4*xcorr(y2.^3,y2)+4*xcorr(y2,y2.^3);
dimt=length(t)-1;
SHG1=shg1(dimt/2:dimt*3/2);
SHG2=shg2(dimt/2:dimt*3/2);
MAC1=mac1(dimt/2:dimt*3/2);
MAC2=mac2(dimt/2:dimt*3/2);
figure(22)
plot(t,MAC1,t,MAC2)

 figure(222)
 plot(t,SHG1/max(SHG1),t,SHG2/max(SHG2))
% 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');xlim([-400 400])