(*  Here I compute N3a and N3b in the 4D case when w' >> w and j =1, ep = 1/10. 
I'll set v0 = 0 and kappa = 1/4, i.e. M = 1 which is what Ray used.
Therefore vH = -4M = -4
I'll also use the scattering parameter F_L that I found an analytic approximation to in the range w = 0.1 to 0.2.  
I'm effectively scaling M out of the equations by defining \bar(w) = w/M and bar(v0) = M v0  so that v0 w  = bar(v0) bar(w).  
  I'll use in this version the variable z = log[w']. 
  *)

(* This is the large w' beta for 2D.  What I need to do is to divide it by F_L in order to turn it into the 4D one for w' >> w.   *)

bHwwp2D = -Exp[-I*(w+wp)*vH]*(Sqrt[w]/(2*Pi*kap*Sqrt[wp]))*Exp[-Pi*w/(2*kap)] * Gamma[I*w/kap]

(* Here is the analytic approximation I have for FL.  Note that I need to check carefully the factors of M that go into this, even though I'm setting M = 1. *)

FL = (-1.1939377639182347`-18.597190132980213` I)+(41.10023317577717` +687.6530825282025` I) x-(592.0381225786749` +14429.639820459033` I) (x^2)+(7111.91712606537` +200125.77678553844` I) (x^3)-(64943.8331654414` +1.9340692528412954`*^6 I) (x^4)+(432088.45699343266` +1.3284352689142372`*^7 I) (x^5)-(2.0587698919373904`*^6+6.4823358495201595`*^7 I) (x^6)+(6.857651819510484`*^6+2.2020032872566673`*^8 I) (x^7)-(1.5194308523785282`*^7+4.959265458817212`*^8 I) (x^8)+(2.0145433966782585`*^7+6.665029117929472`*^8 I) (x^9)-(1.2103820721264469`*^7+4.049821567364419`*^8 I) (x^10);

FL = FL/.{x->w}

(* It works better if I put the a factor of Exp[kap*z/2] here rather than inserting
a factor of Exp[kap*z] directly into the integral over z.  This comes from dwp = kap*z dz.  Also I'll omit the factor of exp(-i vH wp) because it just cancels out
 when you compute Njn.  *)

bHwwpz2D = -Exp[-I*w*vH]*(Sqrt[w]*Exp[-kap*z/2]/(2*Pi*kap^(3/2)))*Exp[-I*w*z]*Exp[-Pi*w/(2*kap)]* Gamma[I*w/kap]

(* It works better if I put the a factor of Exp[z/2] here rather than inserting
a factor of Exp[z] directly into the integral over z.  This comes from dwp = z dz.  I've also inserted the factor of 1/FL  *)

betajnzint2D = (1/Sqrt[ep])*Exp[2*Pi*I*w*n/ep]*bHwwpz2D*kap*Exp[kap*z/2]

betajnzint = (1/Sqrt[ep])*Exp[2*Pi*I*w*n/ep]*bHwwpz2D*kap*Exp[kap*z/2]/FL

BH2D = - Exp[-I*w*vH]*Sqrt[w]*Exp[-Pi*w/(2*kap)]*Gamma[I*w/kap]/(2*Pi*kap)

BH2Ds =  - Exp[I*w*vH]*Sqrt[w]*Exp[-Pi*w/(2*kap)]*Gamma[-I*w/kap]/(2*Pi*kap)

NjnH2Dint = BH2D*BH2Ds

NjnH2D = N[(2*Pi*(1/4)/(1/10))*NIntegrate[NjnH2Dint/.{kap->1/4},{w,1/10,2/10},WorkingPrecision ->15],15]

NjnH2Dalt = N[(2*Pi*(1/4)/(1/10))*NIntegrate[Abs[BH2D/.{vH->-4,kap->1/4}]^2,{w,1/10,2/10},WorkingPrecision ->15],15]

(* This is much much better!  But it might be a useful check to compute the 
4D case both ways.  But for now maybe I'll just do it the easy way! *)

BH4D = BH2D/FL

N3a4D = N[(2*Pi*(1/4)/(1/10))*NIntegrate[Abs[BH4D/.{vH->-4,kap->1/4}]^2,{w,1/10,2/10},WorkingPrecision ->15],15]

(* OK now let's compute N3b by integrating from -Infinity to Lambda for Lambda = 10.  Here is the cutoff in terms of Z  *)

Z = kap*Log[Lam/kap]

ZLam10 =N[Z/.{kap->1/4,Lam->10},15]

(* Recall N3b has a minus sign in front of it. *)

T3ba = Table[{i,-NIntegrate[Abs[NIntegrate[betajnzint/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,-Infinity,-10^5},WorkingPrecision ->10]},{i,1,10}]

T3bb = Table[{i,-NIntegrate[Abs[NIntegrate[betajnzint/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,-10^5,-10^4},WorkingPrecision ->10]},{i,1,10}]

T3bc = Table[{i,-NIntegrate[Abs[NIntegrate[betajnzint/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,-10^4,-10^3},WorkingPrecision ->10]},{i,1,10}]

T3bd = Table[{i,-NIntegrate[Abs[NIntegrate[betajnzint/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,-10^3,-10^2},WorkingPrecision ->10]},{i,1,10}]

T3be = Table[{i,-NIntegrate[Abs[NIntegrate[betajnzint/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,-10^2,-10},WorkingPrecision ->10]},{i,1,10}]

T3bf = Table[{i,-NIntegrate[Abs[NIntegrate[betajnzint/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,-10,0},WorkingPrecision ->10]},{i,1,10}]

T3bg= Table[{i,-NIntegrate[Abs[NIntegrate[betajnzint/.{n->i,ep->1/10,z->zz,kap->1/4,vH->-4},{w,1/10,2/10}]]^2,{zz,0,ZLam10},WorkingPrecision ->10]},{i,1,10}]

TN3b = Table[{i,T3ba[[i,2]]+T3bb[[i,2]]+T3bc[[i,2]]+T3bd[[i,2]]+T3be[[i,2]]+T3bf[[i,2]]+T3bg[[i,2]]},{i,10}]

Tn1N3b = Table[{i,i*TN3b[[i,2]]},{i,10}]

nNjn14D = -((1/4)/(2*Pi))*(Abs[BH4D/.{vH->-4,kap->1/4,w->1/5}]^2+ Abs[BH4D/.{vH->-4,kap->1/4,w->1/10}]^2)

p1 = Plot[nNjn14D,{n,0,11},PlotRange->{{0,11},{-0.00105,-0.00085}}]

p2 = ListPlot[Tn1N3b,PlotRange->{{0,11},{-0.00105,-0.00085}}]

Show[p1,p2]

Export["C:\\Users\\anderson\\Research\\BH_evaporation\\Particle_Production\\2019\\N3b4D.dat",TN3b]