; A demonstrative IDL program to read coefficients and contract rate constants for N2 vibrational excitation/deexcitation 

pro readn2vibcoef

; Impact vibrational exciation led by thermal electrons
 restore, 'impcoef_Tefit.sav'
 
 Ev=[0., 0.2886, 0.5737, 0.8553, 1.1335, 1.4082, 1.6794, 1.947, 2.2111, 2.4717, 2.7287, 2.9821, 3.232, 3.4782, 3.7208, 3.9598, $
   4.1951, 4.4268, 4.6547, 4.8789, 5.0993]*1.6022d4/1.3806 

 Te=4000.d
 Te<=10000.d   ; current polynomial fit is limited to 10000 K
 
 Thele_rt=replicate(0.d,21,21)
; Dimiss_Flag marks transitions with too low rate coefficients and not considered in the model
 for i=0,20 do begin
   for j=i+1,20 do begin
       if Te gt 1700. then begin
          if Dismiss_flag_A[i,j] eq 1 then continue
          t1=alog10(Te/300.d)
          for m=0,6 do Thele_rt[i,j]+=impcoef_Te_A[i,j,m]*(t1^m)
          Thele_rt[i,j]=(10^Thele_rt[i,j])
       endif else if Te ge 500. then begin
         if Dismiss_flag_B[i,j] eq 1 then continue
         t1=Te/300.d
         for m=0,4 do Thele_rt[i,j]=impcoef_Te_B[i,j,m]*(t1^m)
         Thele_rt[i,j]=(10^Thele_rt[i,j])
       endif else if Te ge 200. then begin
          if Dismiss_flag_C[i,j] eq 1 then continue
          t1=Te/300.d
          for m=0, 3 do Thele_rt[i,j]=impcoef_Te_C[i,j,m]*(t1^m)
          Thele_rt[i,j]=(10^Thele_rt[i,j])
       endif

   endfor
 endfor
    
 for i=1, 20 do begin
   for j=0, i-1 do begin
       if Thele_rt[j,i] gt 0. then begin
          Thele_rt[i,j]=Thele_rt[j,i]*exp((Ev[i]-Ev[j])/Te)
       endif else Thele_rt[i,j]=0.d
   endfor
 endfor


; VT exciation led by collision aming N2
 restore, 'VTcoef_n2.sav'

 Tn=1000.d     ; Neutral Temperature
 Tn<=2000.d    ; current polynomial fit is limited to 2000K
 
 VT_N2_rt=replicate(0.d,15,15)
 
 for i=0, 14 do begin
   for j=0, 14 do begin
     if i eq j then continue
     if Tn gt 400. then begin
      if Dismiss_flag_A[i,j] eq 1 then continue
      for n=0, 8 do begin
        VT_N2_rt[i,j]+= VTcoef_N2_A[i,j,n]*(alog10(Tn/300.)^n)
      endfor
      VT_N2_rt[i,j]=10^VT_N2_rt[i,j]
     endif else if Tn ge 200. then begin
      if Dismiss_flag_B[i,j] eq 1 then continue
      for n=0, 3 do begin
        VT_N2_rt[i,j]+= VTcoef_N2_B[i,j,n]*(alog10(Tn/300.)^n)
      endfor
      VT_N2_rt[i,j]=10^VT_N2_rt[i,j]
     endif
   endfor
 endfor
 
 
; VT exciation led by collision with ions
  
 restore, 'VTcoef_i.sav'

 Vdi=4000.d   ; ion flow drift
 Vid <=5000.d ; Current quantum caluclation is limited to 5km/s
 
 Ti=10000.d   ; ion temperature
 bt=1.8046d-3*Vdi*Vdi/Ti   ; see defination in Liang & Donovan (2024)
 bt>=1.       ; current polynomial fit is limited to 1<bt<3
 bt<=3.
 
 VT_i_rt=replicate(0.d,15,15)
 for i=0, 14 do begin
   for j=0, 14 do begin
     if i eq j then continue
     if vdi gt 1750. then begin
       if Dismiss_flag_A[i,j] eq 1 then continue
       for n=0, 7 do begin
         for m=0, 5 do begin
           VT_i_rt[i,j]+= VTcoef_i_A[i,j,n,m]*(alog10(vdi/1000.)^n)*(alog10(bt)^m)
         endfor
       endfor
       VT_i_rt[i,j]=10^VT_i_rt[i,j]
     endif else if vdi ge 1000. then begin
       if Dismiss_flag_B[i,j] eq 1 then continue
       for n=0, 1 do begin
         for m=0, 2 do begin
           VT_i_rt[i,j]+= VTcoef_i_B[i,j,n,m]*(alog10(vdi/1000.)^n)*(alog10(bt)^m)
         endfor
       endfor
       VT_i_rt[i,j]=10^VT_i_rt[i,j]
     endif
 endfor
 endfor

 
; Display the caculated rate constants fro vibrational transitions 
 for i=0, 14 do begin
 for j=0, 14 do begin
  if i eq j then continue
   print, i,j,Thele_rt[i,j],VT_N2_rt[i,j],VT_i_rt[i,j]
 endfor
 endfor

end