%Fundamental wavelength
lambda1=780e-9;

%Crystal thickness
L=6e-3;

%Crystal cut-angle
thetadeg=65;

%Crystal material, valid are 'BBO' and 'Calcite'
ctype='BBO';


lambda2=lambda1/2;
c=3e8;
theta=thetadeg*pi/180;

%Sellmeier
if strcmp(ctype, 'BBO')
    %Coefficients for BBO
    Ao=2.7405;
    Bo=0.0184e-12;
    Co=0.0179e-12;
    Do=0.0155e12;
    Ae=2.3730;
    Be=0.0128e-12;
    Ce=0.0156e-12;
    De=0.0044e12;
end

if strcmp(ctype, 'Calcite')
    %Coefficients for Calcite
    Ao=2.69705;
    Bo=0.0192064e-12;
    Co=0.01820e-12;
    Do=0.0151624e12;
    Ae=2.18438;
    Be=0.0087309e-12;
    Ce=0.01018e-12;
    De=0.0024411e12;
end


% 2:nd harmonic is the extraordinary wave
epso=Ao + Bo/(lambda2^2-Co) - Do*lambda2^2;
n2o=sqrt(epso);
epse=Ae + Be/(lambda2^2-Ce) - De*lambda2^2;
n2e=sqrt(epse);

dn2odlambda2 = (1/n2o) * (-Bo*lambda2/(lambda2^2-Co)^2 - Do*lambda2);
dn2edlambda2 = (1/n2e) * (-Be*lambda2/(lambda2^2-Ce)^2 - De*lambda2);


% fundamental is the ordinary wave
eps1=Ao + Bo/(lambda1^2-Co) - Do*lambda1^2; 
n1=sqrt(eps1);
dn1dlambda1 = (1/n1) * (-Bo*lambda1/(lambda1^2-Co)^2 - Do*lambda1);

deltaeps = epse-epso;



%vector of tuning angles
fi=linspace(-pi/2,pi/2,600); 



D=sqrt(epso*epse*(epso+deltaeps.*(cos(theta).^2)-sin(fi).^2)); 

Q2=epso-sin(fi).^2;

Qe=(D-0.5*deltaeps*sin(2*theta).*sin(fi)) ./ (epso+deltaeps*cos(theta).^2);

AxAz= (Qe.*sin(theta).*sin(fi)-(epso-Qe.^2).*cos(theta)) ./ (Q2.*sin(theta)-Qe.*sin(fi).*cos(theta));

%Extraordinary wave propagation angle
alfa=atan(AxAz);

%Extraordinary wave k-vector angle
alfa_k = atan(sin(fi)./ Qe);

%Ordinary wave propagation angle
beta=asin(sin(fi)/n1);


% Refractive index in the direction of the k-vector
n_alfa_k=1./sqrt(cos(alfa_k+theta).^2./epso+sin(alfa_k+theta).^2/epse);
dn_alfa_k = n_alfa_k.^3 .* (dn2odlambda2.*(cos(alfa_k+theta).^2 ./ n2o.^3) + dn2edlambda2.*(sin(alfa_k+theta).^2 ./ n2e.^3));


%Path for the ordinary wave
L1_kristall = L./cos(beta);
L1_luft = L*(tan(alfa)-tan(beta)).*sin(fi);



%Path for extraordinary wave projected on the k-vector
L2= L./cos(alfa).* cos(alfa-alfa_k);


%Group velocities for fundamental and harmonic
Vg1=c/n1 * (1 + lambda1/n1 .* dn1dlambda1);
Vg2=c./n_alfa_k .* (1 + lambda2./n_alfa_k .* dn_alfa_k);

%Timedelay
tdelay=( L1_kristall./Vg1 + L1_luft./c - L2./Vg2)*1e15;

figure(20)
l1=plot(180/pi*fi,tdelay)
xlabel('Angle of incidence / \circ')
ylabel('\Deltat / fs')
title(['Time delay for ' num2str(L*1000) 'mm ' ctype ', \theta = ' num2str(thetadeg) ' \circ'])
set(l1,'LineWidth',1.5);
