IDL program for calculating Req, EB, Ee in radio TDE

  pro reqmn ;; By G. Mou,  EB+Ee minimal energy method, excluding Ep

;; calculating Req and Eeq in radio TDE

   epse=0.1 & epsb=0.1 & ep=11.0/6.0*epsb/epse

   bet=0.02 ;; 0.22. (beta) SHK velocity. Initially set to be 0, but after determining Vshk, reset it and re-run the procedure.

   Gm =1.0/sqrt(1.0-bet^2.0) ;; 

   fV =0.36 & fA=1.0

;; Source distance setup:

   dL =0.063 & z =0.0457  ;; dL is dL,28. for 2018hyz Cendes 2022.

;   dL =0.071 & z =0.051   ;; dL is dL,28. for 2019dsg Cendes 2021.

;   dL =0.030 & z =0.022  ;; dL is dL,28. for 2019azh Goodwin2022

;   dL=0.24   & z =0.151  ;; dL is dL,28. for 2020opy Goodwin2023.

   zp1=z+1.0

   fourc=4.0; 4.0 if considering the "isotropic equivalent number of electrons". Otherwise: 1.0

   print,"   "

   print,"Standard (EB+Ee) minimal energy method. Theory: Barniol Duran2013."

   print,"Code written by G.Mou on 2023.01.06. "

   print,"Note that this method is rough, which causes several times deviation."

   print,"The result also depends on setups of epsilon_e,B."

   print,"   "

   print,format='("Note: z & dL,28=",f6.4,f7.4,", eps_e,B and ep=",3f6.3)',z,dL,epse,epsb,ep ;

   print,format='("Note: fV & fA=",2f6.2,"  Gm=",f5.2)',fV,fA,Gm

;   p=2.35 & Fp=4.98 & lgvp=9.18 ;; for 2018hyz

   read,"Enter p (p>0):",p

   read,"Enter Fp,mJy: ",Fp

   read,"Enter Lg(Vp): ",lgvp

   vp =10.0^(lgvp-10.0); vp,10

   p13=2.0*p+13.0 ;; 2p+13.

   chi=(p-2.)/(p-1.)*epse*1836. ;; \chi_e

   ksi=1.0 ;; Note: We don't need to include ksi for min(EB+Ee) method, so set it to 1.

   gmm=chi*(Gm-1.0) ;; gmm is \gamma_m, see the paragraph above eq 27 in B.Duran2013.

   if(gmm LT 2.) then gmm=2. ;; \gamma_m, this value should be no less than 2.0.

;;;

   print,format='("CHK: p & vp/GHz & gm_m=",f6.2,f7.3,f7.2)',p,vp,gmm ;

;;; When (EB+Ee) is minimal, the radius is Req: 

   Req=1.0e17*(21.8*525^(p-1.0))^(1.0/p13) * gmm^((2.-p)/p13) * (Fp*dL*dL)^((p+6.)/p13)/vp *  $

      zp1^((-3.0*p-19.0)/p13)* fA^((-5.-p)/p13)*fV^(-1./p13)*(fourc*ksi)^(1.0/p13)* 1.0 *Gm^((p+8.)/p13)

;;; Eeq in E40 erg

   Eeq=1.3e8*21.8^((-2.*p-2.)/p13)* 525.^((11.*p-11.)/p13) * gmm^((22.-11.*p)/p13) * (Fp*dL*dL)^((3.*p+14.)/p13)/vp* $

      zp1^((-5.*p-27.)/p13) * fA^((-3.*p-3.)/p13)* fV^((2.*p+2.)/p13)*(fourc*ksi)^(11./p13) * Gm^((16.-5.*p)/p13)

;;; equipartation method outputs:

   print,"  "

   print,format='("==> Req & LgReq=",g9.3,f7.2," ;   Eeq/1e40 & LgEeq=",g9.3,f7.2)',Req,alog10(Req),Eeq,alog10(Eeq)+40.

;;;

   rfct=ep^(1./17.)  ;; R(real)=rfct * Req

   print,format='("rfct = ep^(1/17) =",f6.4)',rfct

   read,"Enter rfct: ",rfct ;; for Manual revision since Req (coefficient 1e17) is rough, while EB is sensitive to Rrl.

   Rrl=Req * rfct ;; this is the real radius. Once the radius is determined, Ee and EB will be done.

;;

   print,"     "

   print,format='("Physics of Emitting Region:")'

   print,format='("==> R(real) & LgRrl=",g10.3,f8.2," <<==")',Rrl,alog10(Rrl)

;; \gamma_e for \nu_p

   gme=525.*Fp*dL*dL/vp^2.* 1.0 /zp1^3. * Gm/fA/(Rrl/1.0e17)^2.

;;; Ee in E40 erg, eq 17 in B.Duran2013 with a correction of fourc*(gmm/gme)^(2-p)

   Ee= fourc* (gmm/gme)^(2.0-p) * 4.4e10*(Fp*dL*dL)^4.0 /vp^7. * 1.0 * zp1^(-11.) * Gm^2./fA^3./(Rrl/1.0e17)^6.

;;; EB in E40 erg, eq 18 in B.Duran2013, unchanged 

   EB= 2.1e6*(Fp*dL*dL)^(-4.)*vp^10.* 1.0 * zp1^14. * fA^4.*fV* (Rrl/1.0e17)^11./Gm^8.

;;; Eshk

   Eshk=Ee/epse

;;;

   print,format='("gme=",f6.1,"  accounting for \nu_p.")',gme

   print,format='(">>  EB/1e40 & LgEB =",g10.3,f8.2)', EB,alog10(EB)+40.0

   print,format='(">>  Ee/1e40 & LgEe =",g10.3,f8.2)', Ee,alog10(Ee)+40.0

   print,format='("==>  E_shk=",g9.3,"+E40 erg.  Lg Eshk=",f6.2," <<==")',Eshk, alog10(Eshk)+40.

   print,format='("CHK: Lg(EB+Ee)=",f6.2,",  Ee/Eshk & EB/Eshk & Ee/EB=",2f6.3,f7.1)',alog10(EB+Ee)+40., Ee/Eshk,EB/Eshk,Ee/EB

;;;

   print,"  "

   read,"Enter shock velocity via \beta=",betav

   Np =Ee/epse*2.0/(betav*3.0e10)^2/1.67e-24 ;; xE40

   Vol=fV*3.14*(Rrl/1.0e10)^3.0/Gm^4.0

   nnp=Np/Vol*1.0e10

   Bmf=sqrt(EB/Vol*1.0e10*25.1)

   print,format='("Np=",g9.3,"+E40  & np=Np/Vol=",g8.3," cm-3, n_CNM=np/4=",g8.3," cm-3")',Np,nnp,nnp/4.0

   print,format='("B/Gauss & Lg(B/Gauss)=",2f8.3)',Bmf,alog10(Bmf)

;;;

   close, /all

;;

 end