Update with optional effects due to area expansion (#4)

* initial changes from PC * clean up calc_c1 * add new keywords; remove redundante statements * clean up some comments * update the readme * add a tests folder for quick checks * update readme * cleanup comments * last docstring cleanup
rice-solar-physics · Sep 21, 2021 · a309ccd · a309ccd
1 parent c1e36b1
commit a309ccd
Show file tree

Hide file tree

Showing 7 changed files with 428 additions and 53 deletions.
diff --git a/README.md b/README.md
@@ -1,5 +1,7 @@
 # EBTEL
+
 [![ascl:1203.007](https://img.shields.io/badge/ascl-1203.007-blue.svg?colorB=262255)](http://ascl.net/1203.007)
+[![DOI](https://zenodo.org/badge/31337216.svg)](https://zenodo.org/badge/latestdoi/31337216)
 
 Enthalpy-Based Thermal Evolution of Loops (EBTEL) model coded in the IDL language by J. A. Klimchuk, S. Patsourakos, and P.J. Cargill.
 
@@ -12,31 +14,41 @@ For more information regarding the EBTEL model see:
 See also [ebtel++](https://github.com/rice-solar-physics/ebtelPlusPlus), a C++ implementation of the EBTEL model that treats the electron and ion populations separately. Note that while the IDL version of EBTEL only solves the single fluid equations, there is little or no difference between the two-fluid and single-fluid models below roughly 5 MK.
 
 ## Word of Caution
+
 In its integration scheme, EBTEL IDL uses a fixed timestep with a default value of 1 second. In order to accurately model a strong, impulsive heating event, it is necessary to resolve the timescales of the relevant physical processes, most notably thermal conduction, which can be on the order of milliseconds for a typical coronal plasma (e.g. [Bradshaw and Cargill, 2013](http://adsabs.harvard.edu/abs/2013ApJ...770...12B)). Thus the EBTEL IDL integration routine with the default timestep may under-resolve portions of the loop's evolution, particularly in the early heating and conductive cooling phases. 
 
-[ebtel++](https://github.com/rice-solar-physics/ebtelPlusPlus) addresses this issue by providing an adaptive timestep routine that ensures the timestep is always sufficiently small compared to the timescales of the relevant physical processes. ebtel++ also provides an additional correction to the _c_<sub>1</sub> parameter during the conduction phase and treats the electrons and the ions separately (see appendices of [Barnes et al., 2016](http://adsabs.harvard.edu/abs/2016ApJ...829...31B) for more details). Both of these corrections are likely to be important in the case of strong, impulsive heating. Thus, for sufficiently short heating timescales and large heating rates, use of ebtel++ is recommended.
+[ebtel++](https://github.com/rice-solar-physics/ebtelPlusPlus) addresses this issue by providing an adaptive timestep routine that ensures the timestep is always sufficiently small compared to the timescales of the relevant physical processes. This correction is likely to be important in the case of strong, impulsive heating. Thus, for sufficiently short heating timescales and large heating rates, use of ebtel++ is recommended.
 
 ## Terms of Use
+
 This code is the authorized version of EBTEL dated January 2013. Updates will be made as and when necessary.This code is distributed as is. Modifications by the user are unsupported. If you believe you have encountered an error in the original (i.e. unaltered) code, create an issue or submit a pull request with the requested changes.
 
 Use of EBTEL should be acknowledged by referencing Papers 1 & 2 as listed above.  
 
 ## Purpose
+
 Compute the evolution of spatially-averaged (along the field) loop quantities using simplified equations.  The instantaneous differential emission measure of the transition region is also computed. This version incorporates all the modifications from Paper 2 and is written in modular form for clarity. DEM parts unchanged except for TR pressure correction (see Paper 2).
 
 ## Inputs
+
 + ttime  = time array (s)
 + heat   = heating rate array (erg cm^-3 s^-1)   (direct heating only; the first element, heat(0), determines the initial static equilibrium)
 + length = loop half length (top of chromosphere to apex) (cm)
 
 ## Optional Keyword Inputs
+
 + classical = set to use the UNsaturated classical heat flux
 + dem_old   = set to use old technique of computing DEM(T) in the trans. reg. (weighted average of demev, demcon, and demeq)
 + flux_nt   = energy flux array for nonthermal electrons impinging the chromosphere (input as a positive quantity; erg cm^-2 s^-1)
 + energy_nt = mean energy of the nonthermal electrons in keV (default is 50 keV)
 + rtv       = set to use Rosner, Tucker, Vaiana radiative loss function (Winebarger form)
++ a_c       = normalised area factor for the corona (average) as described in literature (see below). (default is 1)
++ a_0       = normalised area factor for the top of the TR as described in literature (see below). (default is 1)
++ a_tr      = normalised area factor for the TR (average) as described in literature (see below). (default is 1)
++ l_fact    = l_cor/length, coronal fraction of loop length recommended value is 0.85 (default is 1)
 
 ## Outputs
+
 + t (t_a) = temperature array corresponding to time (avg. over coronal section of loop / apex)
 + n (n_a) = electron number density array (cm^-3) (coronal avg. / apex)
 + p (p_a) = pressure array (dyn cm^-2) (coronal avg. /apex)
@@ -51,6 +63,7 @@ Compute the evolution of spatially-averaged (along the field) loop quantities us
 + rad_cor =  coronal radiative loss
 
 ## Correspondence with Variables in ApJ Articles
+
 + Paper 1
  + r1 = c_3
  + r2 = c_2
@@ -61,6 +74,7 @@ Compute the evolution of spatially-averaged (along the field) loop quantities us
  + dem_eq = DEM_se
 
 ## Usage
+
 + To include the transition region DEM:
   + `IDL> ebtel2, ttime, heat, t, n, p, v, dem_tr, dem_cor, logtdem`
 + To exclude the transition region DEM (faster):
@@ -71,6 +85,7 @@ Compute the evolution of spatially-averaged (along the field) loop quantities us
   + `IDL> ebtel2, ttime, heat, t, n, p, v, dem_tr, dem_cor, logtdem, f_ratio, rad_ratio`
 
 ### Example Run
+
 + Define the necessary parameters
     + `IDL> time = findgen(10000)` 										Define time array
     + `IDL> heat = fltarr(10000)` 										Define corresponding heating array
@@ -109,6 +124,7 @@ Compute the evolution of spatially-averaged (along the field) loop quantities us
     + `IDL> plot, logtdem, alog10(dem60), xtit='log T (K)',ytit='log DEM (cm!U-5!N K!U-1!N)', tit='1000-1059 s Integration', xran=[5.,7.], /ynoz`
 
 ## Intensities
+
 For observations in which temperature response function, G(T), has units of DN s^-1 pix^-1 cm^5 and the loop diameter, d, is larger than the pixel dimension, l_pix:
 
 + I_cor_perp = d/(2L) * Int{G(T) * dem_cor(T) * dT}
@@ -118,15 +134,19 @@ For observations in which temperature response function, G(T), has units of DN s
 for lines-of-sight perpendicular and parallel to the loop axis. I_tr_perp assumes that the transition region is thinner than l_pix. See `intensity_ebtel.pro` for more information.
 
 ## Miscellaneous Comments
+
 + A 1 sec time step is generally adequate, but in applications where exceptionally strong conductive cooling is expected (e.g., intense short duration heating events, especially in short loops), a shorter time step may be necessary. If there is a question, users should compare runs with different time steps and verify that there are no significant differences in the results.
 + Runs much more quickly if the transition region DEM is not computed.
 + Speed can be increased by increasing the minimum DEM temperature from 10^4 to, say, 10^5 K or by decreasing the maximum DEM temperature from 10^8.5 to, say, 10^7.5 (search on 450 and 451).
 + The equilibrium base heat flux coefficient of 2/7 is appropriate for uniform heating; a coefficient of 4/7 is more appropriate for apex heating.
 + To have equal amounts of thermal and nonthermal heating:  flux_nt = heat*length.
 + It is desirable to have a low-level background heating during the cooling phase so that the coronal temperature does not drop below values at which the corona DEM is invalid.
 + v = (c_2/c_3)(t_tr/t)v_0 = (r2/r1)(t_tr/t)(v/r4) at temperature t_tr in the transition region, where t is the average coronal temperature.
++ The non-uniform area version assumes A_c >= A_0 >= A_TR.  There are limits on how large A_c can be. Too large area factors lead to violation of EBTEL assumptions, in particular that TR thickness is small compared to the loop length. An initial recommendation is A_c < 5. See Cargill et al, MNRAS, 2021 for more information.
 
 
 ## CHANGELOG
+
++ September 2021,  Modified to include non-uniform loop area (see above for caveats) 
 + May 2012. PC version. Modular form.
 + 2013 Jan 15, JAK, Fixed a bug in the compution of the radiation function at 10^4 K; important for computing radiation losses based on dem_tr; ge vs. gt in computing rad;  lt vs. le in computing rad_dem
diff --git a/calc_c1.pro b/calc_c1.pro
@@ -1,13 +1,12 @@
 pro calc_c1, temp, den, length, rad, c1
     common params, k_b, mp, kappa_0
+    common area, a0_ac, atr_ac, l_c, l_tr, l_star
 
     ; Calculate scale height
     calc_lambda, temp, sc
-    ; r3_rad_0: Radiative phase value, no gravity
-    ; r3_eqm_0: Value in eqm with no gravity, -2/3 loss power law.
-    ;
-    r3_rad_0 = 0.6
-    r3_eqm_0 = 2.
+
+    r3_rad_0 = 0.6  ; Radiative phase value, no gravity
+    r3_eqm_0 = 2.  ; Value in eqm with no gravity, -2/3 loss power law.
     l_fact_eq  = 5.
     l_fact_rad = 5.
 
@@ -17,8 +16,8 @@ pro calc_c1, temp, den, length, rad, c1
 
     ; Adjust values for gravity
     ;
-    r3_eqm_g =   r3_eqm_0*exp(2*2*sin(3.14159/l_fact_eq)*length/3.14159/sc)
-    r3_radn_g = r3_rad_0*exp(2*2*sin(3.14159/l_fact_rad)*length/3.14159/sc)
+    r3_eqm_g =   r3_eqm_0*exp(2*2*sin(3.14159/l_fact_eq)*l_c/3.14159/sc)
+    r3_radn_g = r3_rad_0*exp(2*2*sin(3.14159/l_fact_rad)*l_c/3.14159/sc)
     ;
     ; Adjust for loss function
     ;
@@ -27,20 +26,21 @@ pro calc_c1, temp, den, length, rad, c1
     ; These lines turn off loss function fix
 ;   r3_eqm = r3_eqm_g
 ;   r3_radn = r3_radn_g
+    ;
     ;Calculate over/under density
-
-    n_eq_2 = kappa_0/3.5/r3_eqm/rad/length/length*(temp/r2)^(7./2.)
+    ;    
+    n_eq_2 = (a0_ac/atr_ac)*kappa_0/3.5/r3_eqm/rad*(l_star/length/l_c^2)*(temp/r2)^(7./2.)/(1.-l_tr/l_c/r3_eqm)
+
     noneq2=den^2./n_eq_2
-
-    ; Hot loops: use eqm C1
+
+    r3_cond = 6.
+    ;     r3_cond = 2.  ; uncomment this to turn off the conductive phase fix
     if (noneq2 lt 1.) then begin
-    r3=r3_eqm
-    endif
-
-    ; Radiative loops: transition from eqm.
-    if (noneq2 ge 1.) then begin
-    r3=(2.*r3_eqm+r3_radn*(noneq2-1.))/(1+noneq2)
-    endif
+        r3=(2.*r3_eqm+r3_cond*(1./noneq2-1.))/(1+1./noneq2)
+    endif else begin
+        ; Radiative loops: transition from eqm.
+        r3=(2.*r3_eqm+r3_radn*(noneq2-1.))/(1+noneq2)
+    endelse
 
     ; Final value
     c1=r3

diff --git a/ebtel2.pro b/ebtel2.pro
@@ -1,7 +1,8 @@
    pro ebtel2, ttime, heat, length, t, n, p, v, ta, na, pa, c11, dem_tr, dem_cor,  $
        logtdem, f_ratio, rad_ratio, cond, rad_cor,$
        classical=classical, dynamic=dynamic, dem_old=dem_old, $
-       flux_nt=flux_nt, energy_nt=energy_nt, rtv=rtv
+       flux_nt=flux_nt, energy_nt=energy_nt, rtv=rtv, $
+       a_c=a_c, a_0=a_0, a_tr=a_tr, l_fact=l_fact
      ;
      ; NAME:  Enthalpy-Based Thermal Evolution of Loops (EBTEL)
      ;
@@ -22,31 +23,39 @@
      ; OPTIONAL KEYWORD INPUTS:
      ;   classical = set to use the UNsaturated classical heat flux
      ;   dynamic   = set to use dynamical r1 and r2 (NOT recommended, especially when T > 10 MK). Now redundant.
-     ;   dem_old   = set to use old technique of computing DEM(T) in the trans. reg.
+     ;   dem_old   = set to use old technique of computing DEM(T) in the transition region
      ;               (weighted average of demev, demcon, and demeq)
      ;   flux_nt   = energy flux array for nonthermal electrons impinging the chromosphere
      ;               (input as a positive quantity; erg cm^-2 s^-1)
      ;   energy_nt = mean energy of the nonthermal electrons in keV (default is 50 keV)
      ;   rtv       = set to use Rosner, Tucker, Vaiana radiative loss function (Winebarger form)
+     ;   a_c       = normalised area factor for the corona (average) as
+     ;               described in literature (see below). (default is 1)
+     ;   a_0       = normalised area factor for the top of the TR as
+     ;               described in literature (see below). (default is 1)
+     ;   a_tr      = normalised area factor for the TR (average) as
+     ;               described in literature (see below). (default is 1)
+     ;   l_fact    = l_cor/length, coronal fraction of loop length recommended value is 0.85
+     ;               (default is 1)
      ;
      ; OUTPUTS:
-     ;   t (t_a) = temperature array corresponding to time (avg. over coronal section of loop / apex)
-     ;   n (n_a) = electron number density array (cm^-3) (coronal avg. / apex)
-     ;   p (p_a) = pressure array (dyn cm^-2) (coronal avg. /apex)
-     ;   v = velocity array (cm s^-1) (r4 * velocity at base of corona)
-     ;   c11 = C1 (or r3 in this code)
-     ;   dem_tr = differential emission measure of transition region, dem(time,T), both legs
-     ;             (dem = n^2 * ds/dT  cm^-5 K^-1)
-     ;             (Note:  dem_tr is not reliable when a nonthermal electron flux is used.)
-     ;   dem_cor = differential emission measure of corona, dem(time,T), both legs
-     ;             (Int{dem_cor+dem_tr dT} gives the total emission measure of a
-     ;             loop strand having a cross sectional area of 1 cm^2)
-     ;   logtdem = logT array corresponding to dem_tr and dem_cor
-     ;   f_ratio = ratio of heat flux to equilibrium heat flux
-     ;             (ratio of heat flux to tr. reg. radiative loss rate)
+     ;   t (t_a)   = temperature array corresponding to time (avg. over coronal section of loop / apex)
+     ;   n (n_a)   = electron number density array (cm^-3) (coronal avg. / apex)
+     ;   p (p_a)   = pressure array (dyn cm^-2) (coronal avg. /apex)
+     ;   v         = velocity array (cm s^-1) (r4 * velocity at base of corona)
+     ;   c11       = C1 (or r3 in this code)
+     ;   dem_tr    = differential emission measure of transition region, dem(time,T), both legs
+     ;               (dem = n^2 * ds/dT  cm^-5 K^-1)
+     ;               (Note:  dem_tr is not reliable when a nonthermal electron flux is used.)
+     ;   dem_cor   = differential emission measure of corona, dem(time,T), both legs
+     ;               (Int{dem_cor+dem_tr dT} gives the total emission measure of a
+     ;               loop strand having a cross sectional area of 1 cm^2)
+     ;   logtdem   = logT array corresponding to dem_tr and dem_cor
+     ;   f_ratio   = ratio of heat flux to equilibrium heat flux
+     ;               (ratio of heat flux to tr. reg. radiative loss rate)
      ;   rad_ratio = ratio of tr. reg. radiative loss rate from dem_tr and from r3*(coronal rate)
-     ;   cond = conductive loss from corona
-     ;   rad_cor =  coronal radiative loss
+     ;   cond      = conductive loss from corona
+     ;   rad_cor   = coronal radiative loss
      ;
      ; CORRESPONDENCE WITH VARIABLES IN ASTROPHYSICAL JOURNAL ARTICLES:
      ;          (Klimchuk et al., 2008, ApJ, 682, 1351; Cargill et al., ApJ 2012)
@@ -100,17 +109,37 @@
      ;      the same time that r1 = 0.7 provides a more accurate radiative cooling.
      ;   v = (c_2/c_3)*(t_tr/t)*v_0 = (r2/r1)*(t_tr/t)*(v/r4) at temperature t_tr in the transition
      ;      region, where t is the average coronal temperature.
+     ;   The non-uniform area version assumes A_c >= A_0 >= A_TR.  There are limits on how large A_c
+     ;      can be. Too large area factors lead to violation of EBTEL assumptions, in particular that TR
+     ;      thickness is small compared to the loop length. An initial recommendation is A_c < 5. See
+     ;      Cargill et al, MNRAS, 2021 for more information.
      ;
 
      ; HISTORY:
+     ; September 2021,  Modified to include non-uniform loop area (see above for caveats) 
      ; May 2012. PC version. Modular form.
      ; See original ebtel.pro for many additional comments.
      ; 2013 Jan 15, JAK, Fixed a bug in the compution of the radiation function at 10^4 K;
      ;      important for computing radiation losses based on dem_tr;
      ;      ge vs. gt in computing rad;  lt vs. le in computing rad_dem
      ; ---------------------------------------------------------------------
      common params, k_b, mp, kappa_0
-
+     common area, a0_ac, atr_ac, l_c, l_tr, l_star
+
+     if not keyword_set(a_0) then a_0 = 1. 
+     if not keyword_set(a_c) then a_c = 1. 
+     if not keyword_set(a_tr) then a_tr = 1. 
+     if not keyword_set(l_fact) then l_fact = 1. 
+
+     ; Area ratios
+     a0_ac=a_0/a_c
+     atr_ac=a_tr/a_c
+
+     ; Modified lengths
+     l_c = l_fact*length
+     l_tr=length - l_c
+     l_star = l_c + l_tr*atr_ac
+
      ntot = n_elements(ttime)
 
      ; Physical constants Can comment out Hydrad lines if needed.
@@ -242,11 +271,13 @@
 
      ; Iterate on TT and r3
 
+      atr_a0 = a_tr/a_0
+
      for i=1,100 do begin
        calc_c1, tt_old, nn, length, rad, r3
-       tt_new = r2*(3.5*r3/(1. + r3)*length*length*q(0)/kappa_0)^(2./7.)
+       tt_new = r2*(3.5*(r3-l_tr/l_c)/(1. + r3*atr_ac)*l_c^2.*q(0)*atr_a0/kappa_0)^(2./7.)
        radloss, rad, tt_new, 1, rtv=rtv
-       nn = (q(0)/((1. + r3)*rad))^0.5
+       nn = (l_star/l_c*q(0)/((1. + atr_ac*r3)*rad))^0.5
        err = tt_new - tt_old
        err_n = nn_old - nn
        if abs(err) lt 1e3 then  begin
@@ -258,7 +289,7 @@
      endfor
 
      tt = tt_old
-     nn = (q(0)/((1. + r3)*rad))^0.5
+     nn = (l_star/l_c*q(0)/((1. + atr_ac*r3)*rad))^0.5
 
      ;   If want to fix out of eqm start, e.g. cooling flare. Section 4.2, Paper 3.
 ;        tt=1.e7*r2
@@ -277,7 +308,7 @@
      v(0) = 0.
      ta(0) = t(0)/r2
      calc_lambda, t(0), sc
-     na(0) = n(0)*r2*exp(-2*length*(1.-sin(3.14159/5.))/3.14159/sc);
+     na(0) = n(0)*r2*exp(-2*l_c*(1.-sin(3.14159/5.))/3.14159/sc)
      pa(0) = 2*k_b*na(0)*ta(0)
 
      print, 'Initial static equilibrium'
@@ -318,7 +349,7 @@
 
        ; Thermal conduction flux at base
 
-       f_cl = c1*(t(i)/r2)^3.5/length
+       f_cl = c1*(t(i)/r2)^3.5/l_c
 
        if keyword_set(classical) then begin
          f = f_cl
@@ -343,16 +374,15 @@
        ; Equilibrium thermal conduction flux at base (-R_tr in ApJ paper)
        f_eq = -r3*n(i)*n(i)*rad*length
 
-       ;      pv = 0.4*(f_eq - f)
-       pv = 0.4*(f_eq - f - flux_nt(i))
-       ;      dn = pv*0.5/(r12*k_b*t(i)*length)*dt
-       dn = (pv*0.5/(r12*k_b*t(i)*length) + j_nt(i)/length)*dt
+       pv = 0.4*(atr_ac*f_eq/r3/length*(l_c^2/l_star)*(r3-l_tr/l_c) - a0_ac*f - flux_nt(i))
+
+       dn = (pv*0.5/(r12*k_b*t(i)*l_c) + j_nt(i)/length)*dt
 
        n(i+1) = n(i) + dn
-
-       ;      dp = 2./3.*(q(i) + (1. + 1./r3)*f_eq/length)*dt
-       dp = 2./3.*(q(i) + (1. + 1./r3)*f_eq/length    $
-         - (1. - 1.5*k_b*t(i)/energy_nt)*flux_nt(i)/length)*dt
+         
+       dp = 2./3.*(q(i) + (l_c/l_star)*(atr_ac + 1./r3)*f_eq/length    $
+            - (1. - 1.5*k_b*t(i)/energy_nt)*flux_nt(i)/length)*dt
+
        p(i+1) = p(i) + dp
 
        t(i+1) = p(i+1)/(n(i+1)*2.*k_b)
@@ -367,7 +397,7 @@
        ; calculate apex quantities
        ;
        ta(i+1) = t(i+1)/r2;
-       na(i+1) = n(i+1)*r2*exp(-2.*length*(1.-sin(3.14159/5.))/3.14159/sc);
+       na(i+1) = n(i+1)*r2*exp(-2.*l_c*(1.-sin(3.14159/5.))/3.14159/sc);
        pa(i+1) = 2*k_b*na(i+1)*ta(i+1)
 
        ; Differential emission measure
@@ -480,7 +510,7 @@
 
        endif
 
-       cond(i)=f
+       cond(i)=f*a0_ac
        rad_cor(i)=f_eq/r3
 
      endfor