header = "/int/users/paulrom/" <> oheader

file = Import[header <> "freezeout.dat"];

xL = Length[file];

bigN=1;

rfunc[x_] = 
  Normal[Fit[file[[All, 1]], {Sqrt[xL^2-x^2],1, x, x^2, x^3, x^4,x^5, x^6,x^7,x^8}, x]];

tfunc[x_] = 
  Normal[Fit[file[[All, 2]], {1, x}, x]];

Print["xL=",xL];

dtdxfunc[x_] = D[tfunc[x], x];

drdxfunc[x_] = D[rfunc[x], x];

utfunc = Interpolation[file[[All, 3]]];

urfunc = Interpolation[file[[All, 4]]];

pirf=Interpolation[file[[All,5]]];

pief = Interpolation[file[[All, 6]]];

v:=urfunc[x]/utfunc[x]

Ihat0[aa_,bb_]:=1/Pi NIntegrate[Cosh[aa Cos[phi]] Cos[bb Sin[phi]],{phi,0,Pi}];

Ihat1[aa_,bb_]:=1/Pi NIntegrate[Cos[phi] Sinh[aa Cos[phi]] Cos[bb Sin[phi]],{phi,0,Pi}];

Ihat2[aa_,bb_]:=1/Pi NIntegrate[Cos[2 phi] Cosh[aa Cos[phi]] Cos[bb Sin[phi]],{phi,0,Pi}];

Ihat3[aa_,bb_]:=1/Pi NIntegrate[Cos[3 phi] Sinh[aa Cos[phi]] Cos[bb Sin[phi]],{phi,0,Pi}];

Ghat0[aa_,bb_]:=NIntegrate[Exp[-aa Cosh[eta]] Cos[bb Sinh[eta]],{eta,0,Infinity}]

Ghat1[aa_,bb_]:=NIntegrate[Cosh[eta] Exp[-aa Cosh[eta]] Cos[bb Sinh[eta]],{eta,0,Infinity}];

Ghat2[aa_,bb_]:=NIntegrate[Cosh[2 eta] Exp[-aa Cosh[eta]] Cos[bb Sinh[eta]],{eta,0,Infinity}];

Ghat3[aa_,bb_]:=NIntegrate[Cosh[3 eta] Exp[-aa Cosh[eta]] Cos[bb Sinh[eta]],{eta,0,Infinity}];


Ehere[k_,qside_,m_]:=Sqrt[k^2+qside^2/4+m^2];

Ehere1[k_,qout_,m_]:=Sqrt[(k+qout/2)^2+m^2];

Ehere2[k_,qout_,m_]:=Sqrt[(k-qout/2)^2+m^2];

idI0int[k_,T_,m_,x_]:=(k*BesselK[0,Ehere[k,0,m]/T utfunc[x]] 
  BesselI[1,k/T urfunc[x]] dtdxfunc[x]
- Ehere[k,0,m] BesselK[1,Ehere[k,0,m]/T utfunc[x]]
		       BesselI[0,k/T urfunc[x]] drdxfunc[x]);

idI1sideint[k_,qside_,T_,m_,x_]:=(k*BesselK[0,Ehere[k,qside,m]/T utfunc[x]] 
  Ihat1[k/T urfunc[x],qside 5.06842 rfunc[x]] dtdxfunc[x]
- Ehere[k,qside,m] BesselK[1,Ehere[k,qside,m]/T utfunc[x]]
				  Ihat0[k/T urfunc[x],qside 5.06842 rfunc[x]] drdxfunc[x]);


idI1outint[k_,qout_,T_,m_,x_]:=(k*BesselK[0,(Ehere1[k,qout,m]+Ehere2[k,qout,m])/2/T utfunc[x]-I (Ehere1[k,qout,m]-Ehere2[k,qout,m]) 5.06842 tfunc[x]] 
BesselI[1,k/T urfunc[x]-I qout 5.06842 rfunc[x]] dtdxfunc[x]
- (Ehere1[k,qout,m]+Ehere2[k,qout,m])/2 BesselK[1,(Ehere1[k,qout,m]+Ehere2[k,qout,m])/2/T utfunc[x]-I (Ehere1[k,qout,m]-Ehere2[k,qout,m]) 5.06842 tfunc[x]]
BesselI[0,k/T urfunc[x]-I qout 5.06842 rfunc[x]] drdxfunc[x]);

idI1longint[k_,qlong_,T_,m_,x_]:=(k*
Ghat0[Ehere[k,qlong,m]/T utfunc[x],tfunc[x]*5.06842 qlong] 
BesselI[1,k/T urfunc[x]] dtdxfunc[x]
- Ehere[k,qlong,m] 
Ghat1[Ehere[k,qlong,m]/T utfunc[x],tfunc[x] 5.06842 qlong]
BesselI[0,k/T urfunc[x]] drdxfunc[x]);

idI0[k_,T_,m_]:=NIntegrate[rfunc[x] tfunc[x] 
			      idI0int[k,T,m,x],{x,1,xL}, MinRecursion->4, MaxRecursion -> 20];

idI1side[k_,qside_,T_,m_]:=NIntegrate[rfunc[x] tfunc[x] 
					 idI1sideint[k,qside,T,m,x],{x,1,xL}, MinRecursion->4, MaxRecursion -> 20];

idI1out[k_,qout_,T_,m_]:=NIntegrate[rfunc[x] tfunc[x] 
					 Re[idI1outint[k,qout,T,m,x]],{x,1,xL}, MinRecursion->4, MaxRecursion -> 20];

idI2out[k_,qout_,T_,m_]:=NIntegrate[rfunc[x] tfunc[x] 
					 Im[idI1outint[k,qout,T,m,x]],{x,1,xL}, MinRecursion->4, MaxRecursion -> 20];

idI1long[k_,qlong_,T_,m_]:=NIntegrate[rfunc[x] tfunc[x] 
					 idI1longint[k,qlong,T,m,x],{x,1,xL}, MinRecursion->4, MaxRecursion -> 20];


idCout[k_,qout_,T_,m_]:=Module[{re,t1,t2,t3,t4},
			     t1=idI1out[k,qout,T,m];
			     t2=idI2out[k,qout,T,m];
			     t3=idI0[k-qout/2,T,m];
			     t4=idI0[k+qout/2,T,m];
			     re=(t1^2+t2^2)/(t3*t4);
			     re];

idCside[k_,qside_,T_,m_]:=Module[{re,t1,t2,t3,t4},
			     t1=idI1side[k,qside,T,m];
			     t2=idI0[Sqrt[k^2+qside^2/4.],T,m];
			     re=(t1^2)/(t2^2);
			      re];

idClong[k_,qlong_,T_,m_]:=Module[{re,t1,t2,t3,t4},
			     t1=idI1long[k,qlong,T,m];
			     t2=idI0[Sqrt[k^2+qlong^2/4.],T,m];
			     re=(t1^2)/(t2^2);
			      re];


A:=Sqrt[p^2+m^2]/T utfunc[x]

B:=p/T urfunc[x]

mt:=Sqrt[p^2+m^2]

v:=urfunc[x]/utfunc[x]

L1:=2 p mt v pirf[
      x] (mt (BesselK[0,A]+BesselK[2,A]) BesselI[1,B] drdxfunc[x]-
        p dtdxfunc[x] BesselK[1,A] (BesselI[2,B]+BesselI[0,B]))

L2:=-p^2 pirf[
      x] (2 mt BesselK[1,A] BesselI[2,B] drdxfunc[x]-
        p dtdxfunc[x] BesselK[0,A] (BesselI[3,B]+BesselI[1,B]))

L3:=-v^2 mt^2 pirf[
      x] (mt/2 (BesselK[3,A]+3BesselK[1,A]) BesselI[0,B] drdxfunc[x]-
        p dtdxfunc[x] (BesselK[0,A]+BesselK[2,A]) BesselI[1,B])

L4:=-v^2 p^2 pirf[
      x] (mt BesselK[1,A](BesselI[0,B]-BesselI[2,B]) drdxfunc[x]-
        p dtdxfunc[x]/2 BesselK[0,A](BesselI[1,B]-BesselI[3,B]))

L5:=mt^2 pief[
      x](mt 2 BesselK[1,A] BesselI[0,B] drdxfunc[x]-
        p dtdxfunc[x] 2 BesselK[0,A] BesselI[1,B])

L6:=p^2 pief[
      x] (mt BesselK[1,A] (BesselI[0,B]-BesselI[2,B]) drdxfunc[x]-
        p dtdxfunc[x]/2 BesselK[0,A](BesselI[1,B]-BesselI[3,B]))

L7:=-mt^2 pief[
      x](mt/2 (BesselK[3,A]+3 BesselK[1,A]) BesselI[0,B] drdxfunc[x]-
        p dtdxfunc[x] (BesselK[0,A]+BesselK[2,A]) BesselI[1,B])

Ifunc[pp_,mm_,TT_,
    xx_]:=((L1+L2+L3+L4+L5+L6+L7))/.x\[Rule]xx/.T\[Rule]TT/.m\[Rule]
        mm/.p\[Rule]pp

dI0[p_, T_, m_] := -
  1/4./T^2 NIntegrate[rfunc[x] tfunc[x]Ifunc[p, m, T, x], {x, 1, xL}, MinRecursion->4,MaxRecursion -> 20];

Aside:=Eside/T utfunc[x];

Bside:=p/T urfunc[x];

aCisde:=5.06842*q*rfunc[x];


sideL1:=2 p Eside v pirf[x] (
Eside (BesselK[0,Aside]+BesselK[2,Aside]) Ihat1[Bside,aCisde] drdxfunc[x]-
p dtdxfunc[x] BesselK[1,Aside] (Ihat2[Bside,aCisde]+Ihat0[Bside,aCisde]));

sideL2:=-p^2 pirf[x] (2 Eside BesselK[1,Aside] Ihat2[Bside,aCisde] drdxfunc[x]-
		      p dtdxfunc[x] BesselK[0,Aside] (Ihat3[Bside,aCisde]+Ihat1[Bside,aCisde]));

sideL3:=-v^2 Eside^2 pirf[x] (Eside/2 (BesselK[3,Aside]+3BesselK[1,Aside]) 
Ihat0[Bside,aCisde] drdxfunc[x]-
	p dtdxfunc[x] (BesselK[0,Aside]+BesselK[2,Aside]) Ihat1[Bside,aCisde]);

sideL4:=-v^2 p^2 pirf[
		      x] (Eside BesselK[1,Aside](Ihat0[Bside,aCisde]-Ihat2[Bside,aCisde]) drdxfunc[x]-
			  p dtdxfunc[x]/2 BesselK[0,Aside](Ihat1[Bside,aCisde]-Ihat3[Bside,aCisde]));

  sideL5:=Eside^2 pief[x](Eside 2 BesselK[1,Aside] Ihat0[Bside,aCisde] drdxfunc[x]-
		       p dtdxfunc[x] 2 BesselK[0,Aside] Ihat1[Bside,aCisde]);

sideL6:=p^2 pief[x] (Eside BesselK[1,Aside] (Ihat0[Bside,aCisde]-Ihat2[Bside,aCisde]) drdxfunc[x]-
		     p dtdxfunc[x]/2 BesselK[0,Aside](Ihat1[Bside,aCisde]-Ihat3[Bside,aCisde]));

sideL7:=-Eside^2 pief[x](Eside/2 (BesselK[3,Aside]+3 BesselK[1,Aside]) Ihat0[Bside,aCisde] 
drdxfunc[x]-
		      p dtdxfunc[x] (BesselK[0,Aside]+BesselK[2,Aside]) Ihat1[Bside,aCisde]);

sideIfunc[kk_,qq_,mm_,TT_,
      xx_]:=((sideL1+sideL2+sideL3+sideL4+sideL5+sideL6+sideL7))/.x->xx/.T->TT/.Eside->Ehere[kk,qq,mm]/.q->qq/.p->kk


Aout:=Eout/T utfunc[x]-I Ediff 5.06842 tfunc[x];

Bout:=p/T urfunc[x]-I q 5.06842 rfunc[x];

outL1:=2 p Eout v pirf[
      x] (Eout (BesselK[0,Aout]+BesselK[2,Aout]) BesselI[1,Bout] drdxfunc[x]-
        p dtdxfunc[x] BesselK[1,Aout] (BesselI[2,Bout]+BesselI[0,Bout]))

outL2:=-p^2 pirf[
      x] (2 Eout BesselK[1,Aout] BesselI[2,Bout] drdxfunc[x]-
        p dtdxfunc[x] BesselK[0,Aout] (BesselI[3,Bout]+BesselI[1,Bout]))

outL3:=-v^2 Eout^2 pirf[
      x] (Eout/2 (BesselK[3,Aout]+3BesselK[1,Aout]) BesselI[0,Bout] drdxfunc[x]-
        p dtdxfunc[x] (BesselK[0,Aout]+BesselK[2,Aout]) BesselI[1,Bout])

outL4:=-v^2 p^2 pirf[
      x] (Eout BesselK[1,Aout](BesselI[0,Bout]-BesselI[2,Bout]) drdxfunc[x]-
        p dtdxfunc[x]/2 BesselK[0,Aout](BesselI[1,Bout]-BesselI[3,Bout]))

outL5:=Eout^2 pief[
      x](Eout 2 BesselK[1,Aout] BesselI[0,Bout] drdxfunc[x]-
        p dtdxfunc[x] 2 BesselK[0,Aout] BesselI[1,Bout])

outL6:=p^2 pief[
      x] (Eout BesselK[1,Aout] (BesselI[0,Bout]-BesselI[2,Bout]) drdxfunc[x]-
        p dtdxfunc[x]/2 BesselK[0,Aout](BesselI[1,Bout]-BesselI[3,Bout]))

outL7:=-Eout^2 pief[
      x](Eout/2 (BesselK[3,Aout]+3 BesselK[1,Aout]) BesselI[0,Bout] drdxfunc[x]-
        p dtdxfunc[x] (BesselK[0,Aout]+BesselK[2,Aout]) BesselI[1,Bout])

  outIfunc[kk_,qq_,mm_,TT_,
	   xx_]:=((outL1+outL2+outL3+outL4+outL5+outL6+outL7))/.x\[Rule]xx/.T\[Rule]TT/.Eout->(Ehere1[kk,qq,mm]+Ehere2[kk,qq,mm])/2./.q->qq/.p->kk/.Ediff->(Ehere1[kk,qq,mm]-Ehere2[kk,qq,mm])


Along:=Elong/T utfunc[x];

Blong:=p/T urfunc[x];

Cistlong:=tfunc[x]*5.06842 q;

longL1:=2 p Elong v pirf[
			 x] (Elong (Ghat0[Along,Cistlong]+Ghat2[Along,Cistlong]) BesselI[1,Blong] drdxfunc[x]-
        p dtdxfunc[x] Ghat1[Along,Cistlong] (BesselI[2,Blong]+BesselI[0,Blong]))

longL2:=-p^2 pirf[
      x] (2 Elong Ghat1[Along,Cistlong] BesselI[2,Blong] drdxfunc[x]-
        p dtdxfunc[x] Ghat0[Along,Cistlong] (BesselI[3,Blong]+BesselI[1,Blong]))

longL3:=-v^2 Elong^2 pirf[
      x] (Elong/2 (Ghat3[Along,Cistlong]+3Ghat1[Along,Cistlong]) BesselI[0,Blong] drdxfunc[x]-
        p dtdxfunc[x] (Ghat0[Along,Cistlong]+Ghat2[Along,Cistlong]) BesselI[1,Blong])

longL4:=-v^2 p^2 pirf[
      x] (Elong Ghat1[Along,Cistlong](BesselI[0,Blong]-BesselI[2,Blong]) drdxfunc[x]-
        p dtdxfunc[x]/2 Ghat0[Along,Cistlong](BesselI[1,Blong]-BesselI[3,Blong]))

longL5:=Elong^2 pief[
      x](Elong 2 Ghat1[Along,Cistlong] BesselI[0,Blong] drdxfunc[x]-
        p dtdxfunc[x] 2 Ghat0[Along,Cistlong] BesselI[1,Blong])

longL6:=p^2 pief[
      x] (Elong Ghat1[Along,Cistlong] (BesselI[0,Blong]-BesselI[2,Blong]) drdxfunc[x]-
        p dtdxfunc[x]/2 Ghat0[Along,Cistlong](BesselI[1,Blong]-BesselI[3,Blong]))

longL7:=-Elong^2 pief[
      x](Elong/2 (Ghat3[Along,Cistlong]+3 Ghat1[Along,Cistlong]) BesselI[0,Blong] drdxfunc[x]-
	 p dtdxfunc[x] (Ghat0[Along,Cistlong]+Ghat2[Along,Cistlong]) BesselI[1,Blong]);

longIfunc[kk_,qq_,mm_,TT_,
    xx_]:=((longL1+longL2+longL3+longL4+longL5+longL6+longL7))/.x\[Rule]xx/.T\[Rule]TT/.Elong->Ehere[kk,qq,mm]/.q->qq/.p->kk


I1side[k_,qside_,T_,m_] := -1/4./T^2 NIntegrate[
						rfunc[x] tfunc[x] sideIfunc[k, qside, m, T, x], {x, 1, xL}, MinRecursion->4,MaxRecursion -> 20];


I1out[k_,qout_,T_,m_]:=-1/4./T^2 NIntegrate[rfunc[x] tfunc[x] 
				  Re[outIfunc[k,qout,m,T,x]],{x,1,xL}, MinRecursion->4, MaxRecursion -> 20];

I2out[k_,qout_,T_,m_]:=-1/4./T^2 NIntegrate[rfunc[x] tfunc[x] 
				  Im[outIfunc[k,qout,m,T,x]],{x,1,xL}, MinRecursion->4, MaxRecursion -> 20];

I1long[k_,qlong_,T_,m_] := -1/4./T^2 NIntegrate[
						rfunc[x] tfunc[x] longIfunc[k, qlong, m, T, x], {x, 1, xL}, MinRecursion->4,MaxRecursion -> 20];



Cside[k_,qside_,T_,m_]:=Module[{re,t1,t2,t3,t4},
			      t1=idI1side[k,qside,T,m];
			      t3=I1side[k,qside,T,m];
			      t2=idI0[Sqrt[k^2+qside^2/4.],T,m];
			      t4=dI0[Sqrt[k^2+qside^2/4.],T,m];
			     re=(t1+t3)^2/(t2+t4)^2;
			      re];

Cout[k_,qout_,T_,m_]:=Module[{re,t1,t2,t3,t4},
			     t1=idI1out[k,qout,T,m]+I1out[k,qout,T,m];
			     t2=idI2out[k,qout,T,m]+I2out[k,qout,T,m];
			     t3=idI0[k-qout/2,T,m]+dI0[k-qout/2,T,m];
			     t4=idI0[k+qout/2,T,m]+dI0[k+qout/2,T,m];
			     re=(t1^2+t2^2)/(t3*t4);
			     re];

Clong[k_,qlong_,T_,m_]:=Module[{re,t1,t2,t3,t4},
			      t1=idI1long[k,qlong,T,m];
			      t3=I1long[k,qlong,T,m];
			      t2=idI0[Sqrt[k^2+qlong^2/4.],T,m];
			      t4=dI0[Sqrt[k^2+qlong^2/4.],T,m];
			     re=(t1+t3)^2/(t2+t4)^2;
			      re];
