\[Lambda]p0 = 0.788;(*Centerwavelength of pump laser*)
\[Lambda]0 = 2*\[Lambda]p0;(*Centerwavelength of laser*)
c = 2.997*10^14;(*Speed of light*)
vp = c/\[Lambda]p0;(* Pump freq.*)
wp = 2*\[Pi]*vp;(* Angular freq.of pump*)
z0 = 20*L;(* Depth of focus for pump*)
A = Sqrt[z0*\[Lambda]p0/\[Pi]];(* Waist size for pump*)
(*dlp=5*10^-2;*)(* bandwidth for pump*)
dlp = 2.2*10^-2;(* bandwidth for pump*)
dvp = (dlp*c)/\[Lambda]p0^2;
dwp = 2*\[Pi]*dvp;
v = c/\[Lambda]0;(* Center freq.of laser*)
w = 2*\[Pi]*v;(* Angular freq.of laser*)
L = 20*10^3;(* Crystal length*)
(*ordinary refractive index is n2 from modified Fan's equation, \
extrordinary refractive index is n3 from Fradkin's equation*)
n2[\[Lambda]_, T_] :=
Sqrt[2.09930 + 0.922683/(1 - 0.0467695/\[Lambda]^2) -
0.0138404 \[Lambda]^2] + (0.26486/\[Lambda]^3 -
0.60629/\[Lambda]^2 + 0.63061/\[Lambda] + 0.62897) (T -
25) 10^(-5) + (1.3470/\[Lambda]^3 - 3.5770/\[Lambda]^2 +
2.2244/\[Lambda] - 0.14445) (T - 25)^2* 10^(-8);
n3[\[Lambda]_, T_] :=
Sqrt[2.12725 + 1.18431/(1 - 0.0514852/\[Lambda]^2) + 0.6603/(
1 - 100.00507/\[Lambda]^2) -
0.00968956 \[Lambda]^2] + (0.41010/\[Lambda]^3 -
0.89603/\[Lambda]^2 + 0.99228/\[Lambda] + 0.99587) (T -
25) 10^(-5) + (3.1481/\[Lambda]^3 - 9.8136/\[Lambda]^2 +
10.459/\[Lambda] - 1.1882) (T - 25)^2* 10^(-8);
\[Beta]2[\[Lambda]_, T_] :=
1/c (n2[\[Lambda],
T] - \[Lambda] (D[n2[x, T], x] /. x -> \[Lambda]));
\[Beta]3[\[Lambda]_, T_] :=
1/c (n3[\[Lambda],
T] - \[Lambda] (D[n3[x, T], x] /. x -> \[Lambda]));
\[Gamma]2[\[Lambda]_,
T_] := -(\[Lambda]^3/(
2 \[Pi] c^2)) (D[n2[x, T], {x, 2}] /. x -> \[Lambda]);
\[Gamma]3[\[Lambda]_,
T_] := -(\[Lambda]^3/(
2 \[Pi] c^2)) (D[n3[x, T], {x, 2}] /. x -> \[Lambda]);
(*for type II, pump beam is on z-axis, signal beam on y-axis and \
idler on z-axis, the period of the pole is estimated as following*)
polset = 46.146;
Kc = 2*\[Pi]/polset;
ls = 0.06;(*Signal wavelength range for evaluation*)
li = 0.06;(*Idler wavelength range for evaluation*)
(*Pump envelope*)
(*\[Alpha][\[Lambda]s_,\[Lambda]i_]:=Exp[-2Log[2]((2 \[Pi] c)/\
\[Lambda]s +(2 \[Pi] c)/\[Lambda]i -wp)^2/ (dwp)^2];*)
\[Alpha][\[Lambda]s_, \[Lambda]i_] :=
Exp[-2 Log[
2] ((1/\[Lambda]s + 1/\[Lambda]i )^-1 - \[Lambda]p0)^2/ (dlp)^2];
(*Phase matching function*)
(*\[CapitalDelta][\[Lambda]s_,\[Lambda]i_]:=(w/c(2*n2[\[Lambda]0/2]-\
n2[\[Lambda]0]-n3[\[Lambda]0])-Kc)+(((2 \[Pi] \
c)/\[Lambda]s-w)(\[Beta]2[\[Lambda]0]-\[Beta]2[\[Lambda]p0])+((2 \
\[Pi] c)/\[Lambda]i-w)(\[Beta]3[\[Lambda]0]-\[Beta]2[\[Lambda]p0]))+(\
1/2((2 \[Pi] c)/\[Lambda]s-w)^2 \[Gamma]2[\[Lambda]0]+1/2((2 \[Pi] \
c)/\[Lambda]i-w)^2 \[Gamma]3[\[Lambda]0]-1/2((2 \[Pi] \
c)/\[Lambda]s+(2 \[Pi] c)/\[Lambda]i-wp)^2 \[Gamma]2[\[Lambda]p0]);*)
Temp = 33;
\[Gamma]s = \[Beta]2[\[Lambda]0, Temp] - \[Beta]2[\[Lambda]p0, Temp];
\[Gamma]i = \[Beta]3[\[Lambda]0, Temp] - \[Beta]2[\[Lambda]p0, Temp];
poling = wp/c*n2[\[Lambda]p0, Temp] - w/c*n2[\[Lambda]0, Temp] -
w/c*n3[\[Lambda]0, Temp] + Kc;
\[CapitalDelta][\[Lambda]s_, \[Lambda]i_] :=
poling + (((2 \[Pi] c)/\[Lambda]s - w) \[Gamma]s + ((
2 \[Pi] c)/\[Lambda]i - w) \[Gamma]i) + (1/
2 ((2 \[Pi] c)/\[Lambda]s - w)^2 \[Gamma]2[\[Lambda]0, Temp] +
1/2 ((2 \[Pi] c)/\[Lambda]i - w)^2 \[Gamma]3[\[Lambda]0, Temp] -
1/2 ((2 \[Pi] c)/\[Lambda]s + (2 \[Pi] c)/\[Lambda]i -
wp)^2 \[Gamma]2[\[Lambda]p0, Temp]);
Phase[\[Lambda]s_, \[Lambda]i_] :=
Sinc[(\[CapitalDelta][\[Lambda]s, \[Lambda]i] L)/2];
(*Joint spetral amplitude*)
AJ[\[Lambda]s_, \[Lambda]i_] := \[Alpha][\[Lambda]s, \[Lambda]i] \
Phase[\[Lambda]s, \[Lambda]i];
(*group delay btw. signal and idler after passing through PPKTP*)
\[CapitalDelta]t =
L*(\[Beta]3[\[Lambda]0, Temp] - \[Beta]2[\[Lambda]0, Temp])
(*The following defines the three parts of the coincidnece \
integration kernel *)
(*The rotating angle before the PBS*)
(*F1[\[Lambda]s_,\[Lambda]i_,\[Tau]_]:=AJ[\[Lambda]s,\[Lambda]i]^2 \
(1+Sin[4 \[CapitalTheta]] (Cos[(2 \[Pi] c)/\[Lambda]i \
(\[CapitalDelta]t+\[Tau])]+Cos[(2 \[Pi] c)/\[Lambda]s \
(\[CapitalDelta]t+\[Tau])])+Sin[4 \[CapitalTheta]]^2 Cos[(2 \[Pi] c)/\
\[Lambda]i (\[CapitalDelta]t+\[Tau])]Cos[(2 \[Pi] c)/\[Lambda]s (\
\[CapitalDelta]t+\[Tau])]);
F2[\[Lambda]s_,\[Lambda]i_,\[Tau]_]:=AJ[\[Lambda]i,\[Lambda]s]^2 \
(1-Sin[4 \[CapitalTheta]] (Cos[(2 \[Pi] c)/\[Lambda]i \
(\[CapitalDelta]t+\[Tau])]+Cos[(2 \[Pi] c)/\[Lambda]s \
(\[CapitalDelta]t+\[Tau])])+Sin[4 \[CapitalTheta]]^2 Cos[(2 \[Pi] c)/\
\[Lambda]i (\[CapitalDelta]t+\[Tau])]Cos[(2 \[Pi] c)/\[Lambda]s (\
\[CapitalDelta]t+\[Tau])]);
F3[\[Lambda]s_,\[Lambda]i_,\[Tau]_]:=2 AJ[\[Lambda]s,\[Lambda]i]AJ[\
\[Lambda]i,\[Lambda]s](Cos[2 \[CapitalTheta]]^4 Cos[((2 \[Pi] c)/\
\[Lambda]s -(2 \[Pi] c)/\[Lambda]i )(2\[CapitalDelta]t+\[Tau])]+Sin[2 \
\[CapitalTheta]]^4 Cos[((2 \[Pi] c)/\[Lambda]s -(2 \[Pi] \
c)/\[Lambda]i )\[Tau]]-Sin[2 \[CapitalTheta]]^2Cos[2 \
\[CapitalTheta]]^2 (Cos[2*(2 \[Pi] c)/\[Lambda]s\[CapitalDelta]t+((2 \
\[Pi] c)/\[Lambda]s +(2 \[Pi] c)/\[Lambda]i )\[Tau]]+Cos[2*(2 \[Pi] \
c)/\[Lambda]i \[CapitalDelta]t+((2 \[Pi] c)/\[Lambda]s +(2 \[Pi] c)/\
\[Lambda]i )\[Tau]]));*)
(*F3[\[Lambda]s_,\[Lambda]i_,\[Tau]_]:=AJ[\[Lambda]s,\[Lambda]i]AJ[\
\[Lambda]i,\[Lambda]s](Cos[2 \[CapitalTheta]]^4 Cos[((2 \[Pi] c)/\
\[Lambda]s -(2 \[Pi] c)/\[Lambda]i )(\[CapitalDelta]t+\[Tau])]+Sin[2 \
\[CapitalTheta]]^4 Cos[((2 \[Pi] c)/\[Lambda]s -(2 \[Pi] \
c)/\[Lambda]i )(\[CapitalDelta]t-\[Tau])]);*)
F3[\[Lambda]s_, \[Lambda]i_, \[Tau]_] :=
AJ[\[Lambda]s, \[Lambda]i] AJ[\[Lambda]i, \[Lambda]s] (Cos[
2 \[CapitalTheta]]^4 Cos[((2 \[Pi] c)/\[Lambda]s - (
2 \[Pi] c)/\[Lambda]i ) (\[CapitalDelta]t/2 + \[Tau])] +
Sin[2 \[CapitalTheta]]^4 Cos[((2 \[Pi] c)/\[Lambda]s - (
2 \[Pi] c)/\[Lambda]i ) (\[CapitalDelta]t/2 - \[Tau])]);
5.906174747012968`*^-12
P1 = NIntegrate[
Abs[\[Alpha][\[Lambda]s, \[Lambda]i] Phase[\[Lambda]i, \
\[Lambda]s]]^2, {\[Lambda]s, \[Lambda]0 - ls, \[Lambda]0 +
ls}, {\[Lambda]i, \[Lambda]0 - li, \[Lambda]0 + li}]
0.000
13070359382
730715`
P3 = NIntegrate[\[Alpha][\[Lambda]s, \[Lambda]i]^2*
Phase[\[Lambda]s, \[Lambda]i]*
Phase[\[Lambda]i, \[Lambda]s], {\[Lambda]s, \[Lambda]0 -
ls, \[Lambda]0 + ls}, {\[Lambda]i, \[Lambda]0 - li, \[Lambda]0 +
li}]
0.0000879552334912599`
(P3/P1)
0.6729366111192876`
(*P2=NIntegrate[F1[\[Lambda]s,\[Lambda]i,-\[CapitalDelta]t],{\[Lambda]\
s,\[Lambda]0-ls,\[Lambda]0+ls},{\[Lambda]i,\[Lambda]0-li,\[Lambda]0+\
li}]+NIntegrate[F2[\[Lambda]s,\[Lambda]i,-\[CapitalDelta]t],{\[Lambda]\
s,\[Lambda]0-ls,\[Lambda]0+ls},{\[Lambda]i,\[Lambda]0-li,\[Lambda]0+\
li}];
Conincidence1[\[Tau]_]:=1/P1*NIntegrate[F1[\[Lambda]s,\[Lambda]i,\
\[Tau]]+F2[\[Lambda]s,\[Lambda]i,\[Tau]]-F3[\[Lambda]s,\[Lambda]i,\
\[Tau]],{\[Lambda]s,\[Lambda]0-ls,\[Lambda]0+ls},{\[Lambda]i,\[Lambda]\
0-li,\[Lambda]0+li}];*)
P[\[Tau]_] :=
1 + (Sin[4 \[CapitalTheta]]^2 /(2 P1))*
NIntegrate[(AJ[\[Lambda]s, \[Lambda]i]^2 +
AJ[\[Lambda]i, \[Lambda]s]^2) (Cos[((2 \[Pi] c)/\[Lambda]s \
+ (2 \[Pi] c)/\[Lambda]i ) \[Tau]] +
Cos[((2 \[Pi] c)/\[Lambda]s - (2 \[Pi] c)/\[Lambda]i ) \
\[Tau]]) +
AJ[\[Lambda]s, \[Lambda]i]^2 Cos[((2 \[Pi] c)/\[Lambda]s - (2 \
\[Pi] c)/\[Lambda]i ) \[CapitalDelta]t]*
Cos[((2 \[Pi] c)/\[Lambda]s + (2 \[Pi] c)/\[Lambda]i ) \
\[Tau]], {\[Lambda]s, \[Lambda]0 - ls, \[Lambda]0 +
ls}, {\[Lambda]i, \[Lambda]0 - li, \[Lambda]0 + li}] - (1/P1)*
NIntegrate[
F3[\[Lambda]s, \[Lambda]i, \[Tau]], {\[Lambda]s, \[Lambda]0 -
ls, \[Lambda]0 + ls}, {\[Lambda]i, \[Lambda]0 -
li, \[Lambda]0 + li}];
M[\[Tau]_] := \!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Tau]\)]\(P[\[Tau]]\)\);
F[\[Tau]_] = (M[\[Tau]])^2*(1/(P[\[Tau]]*（1 - P[\[Tau]]）));
Plot[F[\[Tau]_], {\[Tau], -1*10^(-12), 1*10^(-12)}, PlotRange -> All]
plot1 = Table[P[\[Tau]], {\[Tau], -1*10^(-12), 1*10^(-12), 1*10^(-13)}]
Save["plot1.dat", plot1]