$HistoryLength=0;

(*Cosmological  Parameters *)
(* CMB Temperature in K *)
TCMB=2.725;
(* Baryon density *)
Ωb=.0458;
(* Cold Dark matter density *)
Ωcdm=.229;
(* Hubble parameter *)
h0=0.702;
(*Helium mass fraction *)
fmHe=0.24;
(* Minimum redshift *)
zmin=200.0;
(* Distortions are calculated below by assuming that all the injected energy goes to radiation (except for baryonic cooling for which electron temperature is evolved).
At z=500 electron temperature differs from radiation temperature by 1% and at z=200 the difference is 10%. 
Below z<200 we cannot asssume that all energy added to the plasma goes to radiation and electron temperature must be evolved
to accurately calculate the y-distortion. See the function dQdzcool[z] below. z<200 is therefore not recommended. *)

(* Maximum redshift *)
zmax=5.0*10^6;
(* Number of neutrinos *)
Nnu=3.046;
(* File with electron fraction with respect to hydrogen (z,ne/nH) as a function of redshift, 
Can be calculated from public codes Recfast,HyRec,CosmoRec with highest redshift  at least before start of helium recombination, z=10000 
and lowest redshift at least 100 *)
xefile="xe.dat";
prefix="test_";
(* Root for template files *)
(*root="intermdediate_templates/";*)
root="intermediate_distortion_templates/";
(*Output file for x,delta n/n*)
outfile1=prefix<>"nxout_nu3.046.dat";

(*output file for ν (GHz), (Subscript[I, ν])(W (m^(-2))(Sr^(-1))(Hz^(-1)))*)
outfile2=prefix<>"Inuout_nu3.046.dat";
resoutfile=prefix<>"resuout_nu3.046.dat";

(*files for error plot output*)
(* Error in total spectrum because of ymin/ymax boundaries for y/mu type distortions *)
error1file=prefix<>"error1.eps";
(* error with respect to numerical results in examples directory *)
error2file=prefix<>"error2.eps";

(*Maxmium and minimum (Subscript[y, γ]) for intermediate type distortions. Above ymax the distortion is assumed to be pure μ and below ymin distortion is pure y*)
ymax=2.0;
ymin=0.01;

(* Check accuracy of final spectrum if checkAccuracy=1 by making ymin=0.001, ymax=10 *)
checkAccuracy=0;
(*Correct for chemical potential suppression at correclow frequencies from double Compton and bremsstrahlung *)
μCorrection=0;
(*Critical frequency to use, μ ∼ μ0*(E^(-xc/x)) or μ0(1.007+3.5*xc)*(((2/π))^(1/2))*(((xc/x))^(1/2))BesselK[0.5*(1-4 *0.7768*xc),xc/x]*)
xc=0.01;
(*Reuse cosmological data: fresh calculation if 0,  otherwise, y parameter and damping scales  are reused from previous runs, saves lot of time*)
reuse=1;
ytozfile="ytoz_reuse.dat";
ztotfile="ztot_reuse.dat";
dampfile="kD_reuse.dat";
tefile="te_reuse.dat";

(*Whether use approximate weights or do the integral to calculate energy release*)
precision=0;
(* Calculate for Silk Damping mode=0, dark matter decay  mode=1, custom energy injection  mode=2 *)
(*If mode is not equal to any of these values, only adiabatic cooling is calculates *)
fisher=1; (*if 0 normal execution, if 1 calculate fisher matrix  *)
fisherfile=prefix<>"fisheroutrunPixiePlanckfid.dat";
(*Add error reduction factor for planck (1.0 for wmap) *)
planck=0;
planckpars=1;
(*Include i-type *)
noitype=0;
(*Experimental parameters for fisher analysis*)
(*Δν Pixie: 15 GHz*)
Δν=15.0;
(*Pixie sensitivity 5*(10^(-26)) W (m^(-2))(Sr^(-1))(Hz^(-1)) *)
ΔI=5.0*(10^(-26));
mode=0;

scalefisher=4;

If[planckpars==1,
(* Hubble parameter *)

h0=0.672;

Ωb=.02225/((h0^2));
(* Cold Dark matter density *)
Ωcdm=.1205/((h0^2));

];


(*Calculate reiduals *)
residual=0;

(*Silk damping parameters for power spectrum *)
(*Spectral index on small scales*)
(* running or not *)
running=0;


If[running==0,
ns=0.96;
(* Amplitude *)
Aζ=1.67*10^-9;
(* pivot point in Mpc*)
k0Mpc=42.0;
dnsdlnk=0.0;,

(*With Running*)
If[planck==1,
(* Hubble parameter *)
ns=0.955;
dnsdlnk=-0.015;
(* Amplitude *)
lnAζ10=3.1;
Aζ=Exp[lnAζ10]*(10^(-10));
k0Mpc=0.05;
,
ns=0.965;
dnsdlnk=0.0;
(* Amplitude *)
Aζ=2.43*10^-9;
(* pivot point in Mpc*)
k0Mpc=0.002;
];
];

(* Dark matter decay parameters *)

(* specify either decay redshift (zORt=0) or lifetime τdecay (zORt=1), (Subscript[n, dm])(t)=(Subscript[n, dm])(initial)*(E^(-t/τdecay)), zdecay is defined by , τdecay=((∫_∞)^zdecay)dt, t is proper time. *)
zORt=0;
(*decay redshift*)
zdecay=10^5;
(*decay time *)
(*τdecay=time[zdecay];*)
(*amplitude of distortion= (((Subscript[f, γ])(Subscript[n, X0])(Subscript[m, X])(c^2))/((Subscript[a, R])(TCMB^4)(1+zdecay)))*)
amplitude=(10^(-7));

(* Custom energy release dQ/dz *)
(* Dark matter annihilation example *)
σv=((3*(10^(-27)))/(Ωcdm (h0^2)));
mdm=10 GeV/(c^2);
ndm[z_]=ρcr*Ωcdm/mdm (((1+z))^3);
fγ=1;
dQdzCustom[z_]= fγ ((mdm (c^2)((ndm[z])^2) σv)/(a (TCMB^4)(((1+z))^5)H[z]));

(*Physical constants *)
(* Compton scattering crossections *)
σT=6.65246*10^(-25);
σ[x_]:=σT*3/4 ((1.+x)/x^3 ((2 x (1+x))/(1+2 x)-Log[1+2 x])+1/(2 x) Log[1+2 x]-(1+3 x)/(1+2 x)^2);
G=6.6742*10^(-8);(*cgs*)
a=7.5658*10^(-15);(*radiation constant cgs*)
c=2.99792458*10^(10);(*cm/s*)
parsec=3.0856776*10^(18);(*cm*)
mparsec=parsec*10^6;
αfs=1/137.036;(* Fine structure constant *)
kb=1.3806505*10^(-16);(*Boltzmann cgs*)
hb=1.05457168*10^(-27);(*Planck/2π cgs*)
h=2*π*hb;
me=9.1093826*10^(-28);(*mass of electron g*)
ev=1.60217646*(10^(-12));
GeV=ev*(10^9);
(*Helium fraction Y=0.24*)
(* Total Matter density *)

Ωm=Ωcdm+Ωb;
H0=h0*100*3.24077649*10^(-20) ;(*s^-1*) 
fHe=fmHe/(3.9715*(1-fmHe));
fnu=Nnu*7.0/8 (4/11)^(4/3);
zeq=(3 (H0 *c)^2 * Ωm)/(8 π G a TCMB^4 (1+fnu))-1;
Ωr=Ωm/(1+zeq);

If[mode==0,zminint=1000.0,zminint=zmin];


H[x_]:=H0*Sqrt[(1+x)^4 * Ωr+Ωm *(1+x)^3+(1-(Ωm+Ωr))];
ρcr=1/(8 π/3*G/H[0]^2);
time[z_]:=NIntegrate[1/((1+x)*H[x]),{x,z,Max[10^4*z,(10^5)]},Method->{"GlobalAdaptive",Method->{"TrapezoidalRule","RombergQuadrature"->True}}];
nH[z_]:=(1+z)^3*Ωb*.76*ρcr/(1.67262*10^(-24)+9.1*10^(-28));
nHe[z_]:=nH[z]*fHe;

xedata=ReadList[xefile,{Real,Real}];
xein=Interpolation[xedata,InterpolationOrder->2];

xezmax=Max[First[First[xedata]],First[Last[xedata]]];
xhighz=xein[xezmax];
xezmin=Min[First[First[xedata]],First[Last[xedata]]];
If[xezmin>101,Print["Warning: Minimum redshift in recombination history must be <=100"];Print["Input file has minimum redshift z=",xezmin]; Exit[];];
xe[z_]=If[z<xezmax,xein[z],xhighz];
ne[z_]=xe[z]*nH[z];
bigA=kb*σT/(me*c);
yy[z_]:=NIntegrate[bigA*ne[x]*TCMB*(1+x)*1/(H[x]*(1+x)),{x,xezmin,z},Method->{"GlobalAdaptive",Method->{"GaussKronrodRule"},"SingularityDepth"->∞},PrecisionGoal->6,MaxRecursion->15];


izmax=3000;
dz=1.0*Log[zmax/zmin]/izmax;
If[reuse==0,
For[i=1,i<=izmax+1,i++,
    zsamp[i]=Exp[Log[zmin]+(i-1)*dz];
    ysamp[i]=Log[yy[zsamp[i]]];
    tsamp[i]=Log[time[zsamp[i]]];
    ];
ytozdata=Table[Flatten[{ysamp[i],Log[zsamp[i]]}],{i,1,izmax+1}];
ztotdata=Table[Flatten[{Log[zsamp[i]],tsamp[i]}],{i,1,izmax+1}];
Export[ytozfile,ytozdata];
Export[ztotfile,ztotdata];,
ytozdata=ReadList[ytozfile,{Real,Real}];
ztotdata=ReadList[ztotfile,{Real,Real}];
];
zdata=Part[Transpose[ztotdata],1];
tdata=Part[Transpose[ztotdata],2];
ttozdata=Transpose[{tdata,zdata}];

ytozint=Interpolation[ytozdata];
ztotint=Interpolation[ztotdata];
ttozint=Interpolation[ttozdata];

ytoz[y_]=Exp[ytozint[Log[y]]];
ztot[z_]=Exp[ztotint[Log[z]]];
ttoz[z_]=Exp[ttozint[Log[z]]];

Print["Solving ODE for Electron temperature"];
coolcomp[z_]:=8 σT a (TCMB (1+z))^4 xe[z]/(3*me c*(1.0+fHe+xe[z]));

If[reuse==0,
solrho=NDSolve[{D[ρ[z],z]==-1/(H[z]*(1+z)) (coolcomp[z]*(1.0-ρ[z]))+ρ[z]/(1+z),ρ[zmax*1.1]==1.0},ρ,{z,zmin,zmax*1.1},Method->{"BDF","MaxDifferenceOrder"->2,"VariableStepCoefficients"->True},AccuracyGoal->15,MaxStepFraction->10^-5,MaxSteps->10^7,InterpolationOrder->All];
Tetbl=Table[Flatten[{Log[zsamp[i]],Log[1.0-ρ[zsamp[i]]/.solrho]}],{i,1,izmax+1}];
Export[tefile,Tetbl];,
Tetbl=ReadList[tefile,{Real,Real}];];

teint=Interpolation[Tetbl];
tgminuste[z_]=Exp[teint[Log[z]]];


(*Calculate Blackbody Optical depth, see  arXiv:1203.2601 for details*)
dcompton[z_]:=1.005*ne[z]*σT c  4*αfs/(3*π)*(((kb *TCMB*(1+z)/(me*c^2)))^2)*25.976;

kcomp[z_]=c ne[z] σT kb TCMB (1+z)/(me c^2);

brehm[z_]:= 2.99*ne[z]*σT c(( αfs*(nH[z]+4*nHe[z]))/(((24 (π^3)))^(1/2)))*((((((me *c^2))/( kb *TCMB *(1+z)))))^(7/2)) *((((h/(me c))))^3);
ztau=1.96*(10^6);
abr=((brehm[ztau])/(((1+ztau))^(5/2)));
adC=((dcompton[ztau])/(((1+ztau))^5));
aC=((kcomp[ztau])/(((1+ztau))^4));
kC=0.7768;
aH=(Ωr (H0^2));

zbr=(((((25 aH)/(4 (kC^2)aC*abr))))^(2/5));

zdC=(((((25 aH)/(4 (kC^2)aC*adC))))^(1/5));
zϵ=((((abr/adC)))^(2/5));

ϵ=((kC 2 abr )/(5 (Sqrt[((adC aH)/aC)])));
zbr22= (((((aH^(1/2))/(2 abr kC *2.958 ))))^2);
zdc22= (((((3 (aH^(1/2)))/( adC kC *2.958 ))))^(1/3));
tau1[z_]=((((((((1+z)/(1+zdC))))^5)+(((((1+z)/(1+zbr))))^(5/2))))^(1/2))+ϵ Log[(((((1+z)/(1+zϵ))))^(5/4))+(((1+(((((1+z)/(1+zϵ))))^(5/2))))^(1/2))];
tau2[z_]=(((((1+z)/(1+zbr22))))^(1/2))+(((((1+z)/(1+zdc22))))^3);

tau[z_]=1.007 tau1[z]+tau2[z];


xmin=10^(-5);
xmax=120.0;
zzmin=Log[xmin];
zzmax=Log[xmax];
nbb[x_]=1.0/(Exp[x]-1.0);

nszsol[x_]=(x Exp[x])/(Exp[x]-1) (x (Exp[x]+1)/(Exp[x]-1)-4);
dIy[ν_]:=inten*(((ν*(10^9)/freq))^3)*nbb[ν*(10^9)/freq]*nszsol[ν*(10^9)/freq];
dIμ[ν_]:=inten*(((ν*(10^9)/freq))^3)*nbb[ν*(10^9)/freq]*musol[ν*(10^9)/freq];
nG[x_]:=(((E^x)*x )/(((-1+(E^x)))^2));
dIt[ν_]:=inten*(((ν*(10^9)/freq))^3)*nG[ν*(10^9)/freq];


If[μCorrection==1,
mucorr[x_]=(1.007+3.5*xc)*(((2/π))^(1/2))*(((xc/x))^(1/2))BesselK[0.5*(1-4 *0.7768*xc),xc/x];,
mucorr[x_]=1;
];
musol[x_]=(Exp[x])/(Exp[x]-1)( (((π^2) )/(18 Zeta[3]))x-mucorr[x]);
dydz[z_]=kb σT ne[z]TCMB (1+z)/(me c H[z] (1+z));

If[mode==0,
(*Silk Damping setup*)
R[z_] := 3/4 *Ωb/(Ωr*(1 + z));
dkddzdamp[z_] := (c*(1 + z))/(6*H[z]*(1 + R[z]) ne[z]*σT) (R[z]^2/(1 + R[z]) + 16/15);
okd2[z_] := NIntegrate[dkddzdamp[zz], {zz, z, ∞}, Method -> {"GlobalAdaptive", Method -> {"GaussKronrodRule"}, "SingularityDepth" -> ∞}, PrecisionGoal -> 6, MaxRecursion -> 30];
If[reuse==0,
kdtbl = Table[Flatten[{z, Log[(1/okd2[Exp[z]])^(1/2)]}], {z, Log[1000.0], Log[zmax], Log[zmax/1000]/500.0}];
Export[dampfile,kdtbl];,
kdtbl=ReadList[dampfile,{Real,Real}];
];
kdint = Interpolation[kdtbl];
kd[z_] := Exp[kdint[Log[z]]];
k0=k0Mpc/mparsec;
If[running==1,
frun[n_,r_,z_]:= NIntegrate[(((k/k0/mparsec))^(n+(1/2)Log[k/k0/mparsec]*r))Exp[-2 (k^2)/(((mparsec*kd[z]))^2)],{k,0,∞},Method->{"GlobalAdaptive",Method->{"GaussKronrodRule"},"SingularityDepth"->∞},PrecisionGoal->5,MaxRecursion->25,MinRecursion->10];
dQdzfisherint[z_,nsnew_,Aζnew_,run_]:=3/4*4.3 Aζnew*k0/mparsec*dkddzdamp[z]*frun[nsnew,run,z];

dQdzfishertable[ns,Aζ,dnsdlnk]=Table[Flatten[{z,Log[dQdzfisherint[Exp[z],ns,Aζ,dnsdlnk]]}],{z,Log[zminint],Log[zmax],Log[zmax/zminint]/1000.0}];
dQdzfisherinter[ns,Aζ,dnsdlnk]=Interpolation[dQdzfishertable[ns,Aζ,dnsdlnk]];
dQdz[z_]=Exp[dQdzfisherinter[ns,Aζ,dnsdlnk][Log[z]]];,
dQdz[z_]=3/4*4.3 Aζ/k0^(ns-1)dkddzdamp[z] 2^(-(3+ns)/2) kd[z]^(ns+1) Gamma[(ns+1)/2];
];

,
If[mode==1,
(* Particle decay setup*)
If[zORt==0,
τdecay=time[zdecay];,
zdecay=ttoz[τdecay];
];
dQdz[z_]=10^(-7)*(1+zdecay)Exp[-ztot[z]/τdecay]/(H[z]*(1+z)^2*τdecay);,
If[mode==2,
(*Custom energy release *)
dQdz[z_]=dQdzCustom[z];,
dQdz[z_]=0;
];
];
];
(* Take into account adiabatic cooling of baryons *)
dQdzcooly[z_]= 4 kb σT  (TCMB (1+z)) ne[z]*tgminuste[z]/(me c*H[z](1+z));
dQdzcool[z_]=3/2 kb TCMB (1+z)*(nH[z]+nHe[z]+ne[z])/(a (TCMB*(1+z))^4)/(1+z);

ystr="nkomp_ad_y5_itr_1_1901_";

y1=0.0;
k=0;
yout[k]=y1;
y2=0.0;
maxflag=1;
minflag=1;
While[y2<=10.0,
      If[y1>1.0,δy=0.01,δy=0.001];
y2=y1+δy;
k=k+1;
yout[k]=y2;
y1=y2;
If[y2≥ymax && maxflag==1,imaxout=k;maxflag=0;];
If[y2≥ymin && minflag==1,iminout=k;minflag=0;];
];
imax=k;
Print["Total Files=",imax];
(*ymin=0.01*)
Print["ymin=",yout[iminout]];
(*ymax=2*)
Print["ymax=",yout[imaxout]];
freq=kb TCMB/h;
inten=(10^(-3))*2 h (freq^3)/(c^2);
z1=ytoz[yout[iminout]];
If[noitype==1,
z2=5*(10^4);,
z2=ytoz[yout[imaxout]];
];

If[fisher==0,
For[i=iminout,i<=imaxout,i++,
(*Print["i=",i];*)
y2=yout[i];
If[i==imaxout,yip1=y2;,yip1=yout[i+1];];
If[i==iminout,yim1=y2;,yim1=yout[i-1];];
dy=(yip1-yim1)/2.0;
(*Calculate weights =energy/4*(10^5)*)
If[precision==0,
yamplitude[i]=(10^5)*0.25*(dQdz[ytoz[y2]]-dQdzcool[ytoz[y2]])/dydz[ytoz[y2]]*dy;,
yamplitude[i]= (10^5)*0.25*NIntegrate[(dQdz[z]-dQdzcool[z]),{z,ytoz[(y2+yim1)/2.0],ytoz[(yip1+y2)/2.0]},Method->{"GlobalAdaptive",Method->{"GaussKronrodRule"},"SingularityDepth"->∞},PrecisionGoal->9,MaxRecursion->20];
];
file=OpenRead[root<>ystr<>ToString[y2]<>".dat"];
data=ReadList[file,{Real,Real}];
If[i==imaxout, xdata=Part[Transpose[data],1];];
nxdata=Part[Transpose[data],2];
If[i==iminout,nx=yamplitude[i]*nxdata;,
nx=nx+yamplitude[i]*nxdata;
];
Close[file];
];
nxfinal=nx;
(*For silk damping y-distortion is only an order of magnitude estimate, cutoff at recombination*)
youtamp= 0.25*NIntegrate[(dQdz[z]-dQdzcooly[z]),{z,zminint,z1},Method->{"GlobalAdaptive",Method->{"GaussKronrodRule"},"SingularityDepth"->∞},PrecisionGoal->9,MaxRecursion->20];
muoutamp= 1.4*NIntegrate[Exp[-tau[z]]*(dQdz[z]-dQdzcool[z]),{z,z2,zmax},Method->{"GlobalAdaptive",Method->{"GaussKronrodRule"},"SingularityDepth"->∞},PrecisionGoal->9,MaxRecursion->20];
youtdata=Table[youtamp*nszsol[Exp[z]],{z,zzmin,zzmax,(zzmax-zzmin)/1000.0}];
muoutdata=Table[muoutamp*musol[Exp[z]],{z,zzmin,zzmax,(zzmax-zzmin)/1000.0}];
(* 150:975 elements of the array select approximately x>(10^(-4)), x<80 to avoid boundaries of the solution in the templates*)
nxoutdata=Transpose[{Part[Exp[xdata],150;;975],Part[nx,150;;975],Part[youtdata,150;;975],Part[muoutdata,150;;975],Part[nx+youtdata+muoutdata,150;;975]}];
Inuoutdata=Transpose[{Part[Exp[xdata]*freq/(10^9),150;;975],Part[inten*((Exp[xdata])^3)*nx*nbb[Exp[xdata]],150;;975],Part[inten*((Exp[xdata])^3)*youtdata*nbb[Exp[xdata]],150;;975],Part[inten*((Exp[xdata])^3)*muoutdata*nbb[Exp[xdata]],150;;975],Part[inten*((Exp[xdata])^3)*(nx+youtdata+muoutdata)*nbb[Exp[xdata]],150;;975]}];
Export[outfile1,nxoutdata];
Export[outfile2,Inuoutdata];
Print["μ-type distortion: μ=",muoutamp];
Print["y-type distortion: y=",youtamp];
youtdatafinal=youtdata;
muoutdatafinal=muoutdata;
];





If[residual==1,
νmin=30+Δν/2.0;
νmax=500;
Inouttable=Transpose[{Part[Exp[xdata]*freq/(10^9),150;;975],Part[inten*((Exp[xdata])^3)*(nx+muoutdata)*nbb[Exp[xdata]],150;;975]}];
Inoutint=Interpolation[Inouttable];
resdata=Table[{ν,Inoutint[ν]},{ν,νmin,νmax,Δν}];
nerror=IntegerPart[(νmax-νmin)/Δν]+1;
Print["No. of points in the fit=",nerror];
errorarray=Flatten[Table[{1/((ΔI)^2)},{i,1,nerror}]];

nmodel[ν_]=y dIy[ν]+μ dIμ[ν]+t dIt[ν];
lm1=NonlinearModelFit[resdata,nmodel[ν],{t,y,μ},ν];
par=lm1["BestFitParameters"];
yfit=y/.par;
μfit=μ/.par;
tfit=t/.par;
Print["Best fit parameters:t=",tfit, " y=",yfit," μ=",μfit];
nfit[ν_]=yfit dIy[ν]+μfit dIμ[ν]+tfit dIt[ν];
resoutdata=Table[{ν,Inoutint[ν],nfit[ν],Inoutint[ν]-nfit[ν]},{ν,νmin,νmax,Δν}];
Export[resoutfile,resoutdata];
];

If[fisher==1,
Print["Calculating Fisher Matrix"];
If[planck==1,
δlnAζ10=lnAζ10/200.0;
δAζ=(10^(-10))Exp[lnAζ10]*δlnAζ10;,
δAζ=Aζ/200.0;
];

δns=ns/200.0;
δdnsdlnk=Max[Abs[dnsdlnk],(10^(-2))]/200.0;
(*ymin=0.01*)
Print["ymin=",yout[iminout]];
(*ymax=2*)
Print["ymax=",yout[imaxout]];
For[j=-1,j≤1,j=j+2,
If[running==1,
dQdzfishertable[ns+j*δns,Aζ,dnsdlnk]=Table[Flatten[{z,Log[dQdzfisherint[Exp[z],ns+j*δns,Aζ,dnsdlnk]]}],{z,Log[zminint],Log[zmax],Log[zmax/zminint]/1000.0}];
dQdzfisherinter[ns+j*δns,Aζ,dnsdlnk]=Interpolation[dQdzfishertable[ns+j*δns,Aζ,dnsdlnk]];
dQdzfisher[z_,ns+j*δns,Aζ,dnsdlnk]=Exp[dQdzfisherinter[ns+j*δns,Aζ,dnsdlnk][Log[z]]];

dQdzfishertable[ns,Aζ+j*δAζ,dnsdlnk]=Table[Flatten[{z,Log[dQdzfisherint[Exp[z],ns,Aζ+j*δAζ,dnsdlnk]]}],{z,Log[zminint],Log[zmax],Log[zmax/zminint]/1000.0}];
dQdzfisherinter[ns,Aζ+j*δAζ,dnsdlnk]=Interpolation[dQdzfishertable[ns,Aζ+j*δAζ,dnsdlnk]];
dQdzfisher[z_,ns,Aζ+j*δAζ,dnsdlnk]=Exp[dQdzfisherinter[ns,Aζ+j*δAζ,dnsdlnk][Log[z]]];

dQdzfishertable[ns,Aζ,dnsdlnk+j*δdnsdlnk]=Table[Flatten[{z,Log[dQdzfisherint[Exp[z],ns,Aζ,dnsdlnk+j*δdnsdlnk]]}],{z,Log[zminint],Log[zmax],Log[zmax/zminint]/1000.0}];
dQdzfisherinter[ns,Aζ,dnsdlnk+j*δdnsdlnk]=Interpolation[dQdzfishertable[ns,Aζ,dnsdlnk+j*δdnsdlnk]];
dQdzfisher[z_,ns,Aζ,dnsdlnk+j*δdnsdlnk]=Exp[dQdzfisherinter[ns,Aζ,dnsdlnk+j*δdnsdlnk][Log[z]]];
,
dQdzfisher[z_,nsnew_,Aζnew_,run_]=3/4*4.3 Aζnew/k0^(nsnew-1)dkddzdamp[z] 2^(-(3+nsnew)/2) kd[z]^(nsnew+1) Gamma[(nsnew+1)/2];
dQdzfisherint[z_,nsnew_,Aζnew_,run_]=3/4*4.3 Aζnew/k0^(nsnew-1)dkddzdamp[z] 2^(-(3+nsnew)/2) kd[z]^(nsnew+1) Gamma[(nsnew+1)/2];
];
For[i=iminout,i<=imaxout,i++,
(*Print["i=",i];*)
y2=yout[i];
If[i==imaxout,yip1=y2;,yip1=yout[i+1];];
If[i==iminout,yim1=y2;,yim1=yout[i-1];];
dy=(yip1-yim1)/2.0;
(*Calculate weights =energy/4*(10^5)*)
If[precision==0,
yamplitudens[i][j]=(10^5)*0.25*(dQdzfisher[ytoz[y2],ns+j*δns,Aζ,dnsdlnk]-dQdzcool[ytoz[y2]])/dydz[ytoz[y2]]*dy;
yamplitudeAζ[i][j]=(10^5)*0.25*(dQdzfisher[ytoz[y2],ns,Aζ+j*δAζ,dnsdlnk]-dQdzcool[ytoz[y2]])/dydz[ytoz[y2]]*dy;
If[running==1,
yamplitudednsdlnk[i][j]=(10^5)*0.25*(dQdzfisher[ytoz[y2],ns,Aζ,dnsdlnk+j*δdnsdlnk]-dQdzcool[ytoz[y2]])/dydz[ytoz[y2]]*dy;
];,
yamplitudens[i][j]= (10^5)*0.25*NIntegrate[(dQdzfisher[z,ns+j*δns,Aζ,dnsdlnk]-dQdzcool[z]),{z,ytoz[(y2+yim1)/2.0],ytoz[(yip1+y2)/2.0]},Method->{"GlobalAdaptive",Method->{"GaussKronrodRule"},"SingularityDepth"->∞},PrecisionGoal->9,MaxRecursion->20];
yamplitudeAζ[i][j]= (10^5)*0.25*NIntegrate[(dQdzfisher[z,ns,Aζ+j*δAζ,dnsdlnk]-dQdzcool[z]),{z,ytoz[(y2+yim1)/2.0],ytoz[(yip1+y2)/2.0]},Method->{"GlobalAdaptive",Method->{"GaussKronrodRule"},"SingularityDepth"->∞},PrecisionGoal->9,MaxRecursion->20];
If[running==1,
yamplitudednsdlnk[i][j]= (10^5)*0.25*NIntegrate[(dQdzfisher[z,ns,Aζ,dnsdlnk+j*δdnsdlnk]-dQdzcool[z]),{z,ytoz[(y2+yim1)/2.0],ytoz[(yip1+y2)/2.0]},Method->{"GlobalAdaptive",Method->{"GaussKronrodRule"},"SingularityDepth"->∞},PrecisionGoal->9,MaxRecursion->20];
];
];
file=OpenRead[root<>ystr<>ToString[y2]<>".dat"];
data=ReadList[file,{Real,Real}];
If[i==imaxout, xdata=Part[Transpose[data],1];];
nxdata=Part[Transpose[data],2];
If[i==iminout,
nxns[j]=yamplitudens[i][j]*nxdata;
nxAζ[j]=yamplitudeAζ[i][j]*nxdata;
If[running==1,
nxdnsdlnk[j]=yamplitudednsdlnk[i][j]*nxdata;];,
nxns[j]=nxns[j]+yamplitudens[i][j]*nxdata;
nxAζ[j]=nxAζ[j]+yamplitudeAζ[i][j]*nxdata;
If[running==1,
nxdnsdlnk[j]=nxdnsdlnk[j]+yamplitudednsdlnk[i][j]*nxdata;
];
];
Close[file];
];
muoutampns[j]= 1.4*NIntegrate[Exp[-tau[z]]*(dQdzfisher[z,ns+j*δns,Aζ,dnsdlnk]-dQdzcool[z]),{z,z2,zmax},Method->{"GlobalAdaptive",Method->{"GaussKronrodRule"},"SingularityDepth"->∞},PrecisionGoal->9,MaxRecursion->20];
muoutampAζ[j]= 1.4*NIntegrate[Exp[-tau[z]]*(dQdzfisher[z,ns,Aζ+j*δAζ,dnsdlnk]-dQdzcool[z]),{z,z2,zmax},Method->{"GlobalAdaptive",Method->{"GaussKronrodRule"},"SingularityDepth"->∞},PrecisionGoal->9,MaxRecursion->20];
If[running==1,
muoutampdnsdlnk[j]= 1.4*NIntegrate[Exp[-tau[z]]*(dQdzfisher[z,ns,Aζ,dnsdlnk+j*δdnsdlnk]-dQdzcool[z]),{z,z2,zmax},Method->{"GlobalAdaptive",Method->{"GaussKronrodRule"},"SingularityDepth"->∞},PrecisionGoal->9,MaxRecursion->20];
];

muoutdatans[j]=Table[muoutampns[j]*musol[Exp[z]],{z,zzmin,zzmax,(zzmax-zzmin)/1000.0}];
muoutdataAζ[j]=Table[muoutampAζ[j]*musol[Exp[z]],{z,zzmin,zzmax,(zzmax-zzmin)/1000.0}];
If[running==1,
muoutdatadnsdlnk[j]=Table[muoutampdnsdlnk[j]*musol[Exp[z]],{z,zzmin,zzmax,(zzmax-zzmin)/1000.0}];
];
If[noitype==1,

Inuoutdatans[j]=inten*((Exp[xdata])^3)*(muoutdatans[j])*nbb[Exp[xdata]];
InuoutdataAζ[j]=inten*((Exp[xdata])^3)*(muoutdataAζ[j])*nbb[Exp[xdata]];
If[running==1,
Inuoutdatadnsdlnk[j]=inten*((Exp[xdata])^3)*(muoutdatadnsdlnk[j])*nbb[Exp[xdata]];
];,

Inuoutdatans[j]=inten*((Exp[xdata])^3)*(nxns[j]+muoutdatans[j])*nbb[Exp[xdata]];
InuoutdataAζ[j]=inten*((Exp[xdata])^3)*(nxAζ[j]+muoutdataAζ[j])*nbb[Exp[xdata]];
If[running==1,
Inuoutdatadnsdlnk[j]=inten*((Exp[xdata])^3)*(nxdnsdlnk[j]+muoutdatadnsdlnk[j])*nbb[Exp[xdata]];
];
];
];

Export[outfile2,Transpose[{Part[Exp[xdata]*freq/(10^9),150;;975],Part[InuoutdataAζ[1],150;;975]}]];
Print["μ-type distortion: μ=",muoutampAζ[1]];

dIns=(Inuoutdatans[1]-Inuoutdatans[-1])/(2*δns);
If[planck==1,
dIAζ=(InuoutdataAζ[1]-InuoutdataAζ[-1])/(2*δlnAζ10);,
dIAζ=(InuoutdataAζ[1]-InuoutdataAζ[-1])/(2*δAζ);
];
Inuoutdata=Transpose[{Part[Exp[xdata]*freq/(10^9),150;;975],Part[dIns,150;;975],Part[dIAζ,150;;975]}];
Export["Fisherderivatives.005.dat",Inuoutdata];
If[running==1,
dIdnsdlnk=(Inuoutdatadnsdlnk[1]-Inuoutdatadnsdlnk[-1])/(2*δdnsdlnk);
Inuoutdata=Transpose[{Part[Exp[xdata]*freq/(10^9),150;;975],Part[dIns,150;;975],Part[dIAζ,150;;975],Part[dIdnsdlnk,150;;975]}];
Export["Fisherderivativesrunning.005.dat",Inuoutdata];

];
dInsint=Interpolation[Transpose[{Part[Exp[xdata]*freq/(10^9),150;;975],Part[dIns,150;;975]}]];
dIAζint=Interpolation[Transpose[{Part[Exp[xdata]*freq/(10^9),150;;975],Part[dIAζ,150;;975]}]];
If[running==1,
dIdnsdlnkint=Interpolation[Transpose[{Part[Exp[xdata]*freq/(10^9),150;;975],Part[dIdnsdlnk,150;;975]}]];
];
];

If[fisher==1,
fisherfile=prefix<>"fisherpixie.dat";
scalefisher=1;
(*Experimental parameters for fisher analysis*)
(*Δν Pixie: 15 GHz*)
Δν=15.0;
(*Pixie sensitivity 5*(10^(-26)) W (m^(-2))(Sr^(-1))(Hz^(-1)) *)
ΔI=5.0*(10^(-26));

ΔIprism1=3.6*(10^(-27));
ΔIprism2=6.5*(10^(-27));
prism=0;

If[running==1,
(* parameters: 1=ns,3=Aζ,2=dnsdlnk,4=t,5=y,6=μ *)
dI[ν_][1]=dInsint[ν];
dI[ν_][3]=(10^(-9))*dIAζint[ν];
dI[ν_][2]=dIdnsdlnkint[ν];
dI[ν_][4]=(10^(-8))*dIt[ν];
dI[ν_][5]=(10^(-8))*dIy[ν];
dI[ν_][6]=(10^(-8))*dIμ[ν];

If[planck==1,
cmbcov={{ 0.53765*10^-04 ,   0.17285*10^-04 ,   0.13729*10^-04},{ 0.17285*10^-04  ,  0.73019*10^-04  , -0.13078*10^-03},{ 0.13729*10^-04 ,  -0.13078*10^-03   , 0.95029*10^-03}};
MatrixForm[cmbcov];,

cmbcov={{0.87258*10^-3,-0.32081*10^-3,-0.17539*10^-2},{-0.32081*10^-3,0.13361*10^-03,0.49953*10^-03},{-0.17539*10^-02,0.49953*10^-03,0.85208*10^-02}};
MatrixForm[cmbcov];
];
If[prism==1,
cmbcov=1/scalefisher*cmbcov;
];
cmbcovnsrun=cmbcov[[1;;2,1;;2]];
MatrixForm[cmbcovnsrun];

cmbfisher=Inverse[cmbcov];
cmbfishernsrun=Inverse[cmbcovnsrun];

cmbFisherMatrix=Table[0,{i,1,6},{j,1,6}];
For[j=1,j≤3,j++,
For[i=1,i≤3,i++,
cmbFisherMatrix[[i,j]]=cmbfisher[[i,j]];
];];
For[j=1,j≤6,j++,
For[i=1,i≤6,i++,
If[prism==1,
Fisher[i,j]=cmbFisherMatrix[[i,j]]+(1/((ΔIprism1)^2))Sum[dI[ν][i]*dI[ν][j],{ν,30+Δν/2.0,180,Δν}]+(1/((ΔIprism2)^2))Sum[dI[ν][i]*dI[ν][j],{ν,180+Δν/2.0,500,Δν}];
,
Fisher[i,j]=cmbFisherMatrix[[i,j]]+(1/((ΔI)^2))Sum[dI[ν][i]*dI[ν][j],{ν,30+Δν/2.0,500,Δν}];
];
];];
Fishermat1=Array[Fisher,{3,3}];
Fishermat2=Array[Fisher,{4,4}];
Fishermat=Array[Fisher,{5,5}];
Fishermatμ=Array[Fisher,{6,6}];
(*Margianlise to get constraints on ns,dnsdlnk*)

A=Fishermat1[[1;;2,1;;2]];
B=Fishermat1[[1;;2,3;;3]];
Cc=Fishermat1[[3;;3,3;;3]];
Fishermarg1=A- B.Inverse[Cc].Transpose[B];



A=Fishermat2[[1;;2,1;;2]];
B=Fishermat2[[1;;2,3;;4]];
Cc=Fishermat2[[3;;4,3;;4]];
Fishermarg2=A-B.Inverse[Cc].Transpose[B];

A=Fishermat[[1;;2,1;;2]];
B=Fishermat[[1;;2,3;;5]];
Cc=Fishermat[[3;;5,3;;5]];
Fishermarg=A-B.Inverse[Cc].Transpose[B];

A=Fishermatμ[[1;;2,1;;2]];
B=Fishermatμ[[1;;2,3;;6]];
Cc=Fishermatμ[[3;;6,3;;6]];
Fishermargμ=A-B.Inverse[Cc].Transpose[B];



Δθ={θx-ns,θy-dnsdlnk};
constrainteq1[θx_,θy_]=Δθ.Fishermarg1.Δθ;
constrainteq2[θx_,θy_]=Δθ.Fishermarg2.Δθ;
constrainteq[θx_,θy_]=Δθ.Fishermarg.Δθ;
constrainteqμ[θx_,θy_]=Δθ.Fishermargμ.Δθ;
constrainteqcmb[θx_,θy_]=Δθ.cmbfishernsrun.Δθ;


outfile=OpenWrite[fisherfile];
If[prism==1,
Export[outfile,{"#Fisher constraints for sensitivity (30-180)=",ΔIprism1,"(180-500)=",ΔIprism2,"bandwidth=",Δν," GHz\n"}];,

Export[outfile,{"#Fisher constraints for sensitivity=",ΔI, "bandwidth=",Δν," GHz\n"}];
];
Export[outfile,FortranForm[Simplify[constrainteq[x,y]]]]
Export[outfile,"\n"]
Export[outfile,"With mu freedom\n"]
Export[outfile,FortranForm[Simplify[constrainteqμ[x,y]]]]
Export[outfile,"\nWith no y,t uncertanity\n"]
Export[outfile,FortranForm[Simplify[constrainteq1[x,y]]]]
Export[outfile,"\nWith no y uncertanity\n"]
Export[outfile,FortranForm[Simplify[constrainteq2[x,y]]]]
Export[outfile,"\ncmb+SPT+ACT\n"]
Export[outfile,FortranForm[Simplify[constrainteqcmb[x,y]]]]
Close[outfile]
,
(* parameters: 1=ns,2=Aζ,3=t,4=y,5=μ *)
dI[ν_][1]=dInsint[ν];
dI[ν_][2]=(10^(-9))*dIAζint[ν];
dI[ν_][3]=(10^(-8))*dIt[ν];
dI[ν_][4]=(10^(-8))*dIy[ν];
dI[ν_][5]=(10^(-8))*dIμ[ν];
For[j=1,j≤5,j++,
For[i=1,i≤5,i++,
If[prism==1,
Fisher[i,j]=(1/((ΔIprism1)^2))Sum[dI[ν][i]*dI[ν][j],{ν,30+Δν/2.0,180,Δν}]+(1/((ΔIprism2)^2))Sum[dI[ν][i]*dI[ν][j],{ν,180+Δν/2.0,500,Δν}];,
Fisher[i,j]=(1/((ΔI)^2))Sum[dI[ν][i]*dI[ν][j],{ν,30+Δν/2.0,500,Δν}];
];
];];
Fishermat1=Array[Fisher,{2,2}];
Fishermat2=Array[Fisher,{3,3}];
Fishermat=Array[Fisher,{4,4}];
Fishermatμ=Array[Fisher,{5,5}];

Fishermarg1=Fishermat1;

A=Fishermat2[[1;;2,1;;2]];
B=Fishermat2[[1;;2,3;;3]];
Cc=Fishermat2[[3;;3,3;;3]];
Fishermarg2=A- B.Inverse[Cc].Transpose[B];



A=Fishermat[[1;;2,1;;2]];
B=Fishermat[[1;;2,3;;4]];
Cc=Fishermat[[3;;4,3;;4]];
Fishermarg=A-B.Inverse[Cc].Transpose[B];

A=Fishermatμ[[1;;2,1;;2]];
B=Fishermatμ[[1;;2,3;;5]];
Cc=Fishermatμ[[3;;5,3;;5]];
Fishermargμ=A-B.Inverse[Cc].Transpose[B];

detectionlimit=Sqrt[1.0/Fishermarg[[2,2]]];

(*
covariance1=Inverse[Fishermat1];
covariance2=Inverse[Fishermat2];
covariance=Inverse[Fishermat];
covarianceμ=Inverse[Fishermatμ];
(*Margianlise to get constraints on ns,Aζ*)
covmarg1=covariance1[[1;;2,1;;2]];
covmarg2=covariance2[[1;;2,1;;2]];
covmarg=covariance[[1;;2,1;;2]];
covmargμ=covarianceμ[[1;;2,1;;2]];
Fishermarg=Inverse[covmarg];
Fishermarg1=Inverse[covmarg1];
Fishermarg2=Inverse[covmarg2];
Fishermargμ=Inverse[covmargμ];

*)

Δθ={θx-ns,θy-Aζ*(10^9)};
constrainteq1[θx_,θy_]=Δθ.Fishermarg1.Δθ;
constrainteq2[θx_,θy_]=Δθ.Fishermarg2.Δθ;
constrainteq[θx_,θy_]=Δθ.Fishermarg.Δθ;
constrainteqμ[θx_,θy_]=Δθ.Fishermargμ.Δθ;
outfile=OpenWrite[fisherfile];
If[prism==1,
Export[outfile,{"#Fisher constraints for sensitivity (30-180)=",ΔIprism1,"(180-500)=",ΔIprism2,"bandwidth=",Δν," GHz\n"}];,
Export[outfile,{"#Fisher constraints for sensitivity=",ΔI, "bandwidth=",Δν," GHz\n"}];
];
Export[outfile,{"1-sigma Detection limit on amplitude for ns=",ns,"is A_zeta=",FortranForm[1.0*(10^(-9))*detectionlimit]}];
Export[outfile,"\n"]
Export[outfile,FortranForm[Simplify[constrainteq[x,y]]]]
Export[outfile,"\n"]
Export[outfile,"With mu freedom\n"]
Export[outfile,FortranForm[Simplify[constrainteqμ[x,y]]]]
Export[outfile,"\nWith no y,t uncertanity\n"]
Export[outfile,FortranForm[Simplify[constrainteq1[x,y]]]]
Export[outfile,"\nWith no y uncertanity\n"]
Export[outfile,FortranForm[Simplify[constrainteq2[x,y]]]]
Close[outfile]
];
];