diff --git a/radynpy/matsplotlib/ContribFn.py b/radynpy/matsplotlib/ContribFn.py
index 70c3231..508bea9 100644
--- a/radynpy/matsplotlib/ContribFn.py
+++ b/radynpy/matsplotlib/ContribFn.py
@@ -68,7 +68,7 @@ def contrib_fn(cdf, kr, tStep=0, yRange=[-0.08, 2.5], vRange=[300.0, -300.0], mu
Print progress information to stdout. (default: False)
returnData : bool, optional
Return a dictionary of the data computed to produce the plots. (default: False)
- opctab : str, optional
+ opctabPath : str, optional
Path to non-standard opctab.dat if needed. (default: False)
'''
with warnings.catch_warnings():
@@ -90,7 +90,6 @@ def contrib_fn(cdf, kr, tStep=0, yRange=[-0.08, 2.5], vRange=[300.0, -300.0], mu
iel = cdf.ielrad[kr] - 1
# Line intensity data
- outMu = cdf.outint[tStep,:,mu,:]
x_ny, tauq_ny = xt_calc(cdf, tStep, iel, kr, withBackgroundOpacity=withBackgroundOpacity, opctabPath=opctabPath)
dtau = np.zeros((nDep, cdf.nq[kr]))
dtau[1:,:] = tauq_ny[1:nDep,:] - tauq_ny[:nDep-1,:]
@@ -324,7 +323,7 @@ def add_image_label(ax, label):
ax[2,0].plot(cdf.vz1[tStep][iwy] * 1e-5, y[iwy], c=vzColor, ls='--')
ax[2,0].plot(x, tau1, c=tauColor)
# Adjust line profile to fill plot
- lineProfile = cdf.outint[tStep, 1:cdf.nq[kr]+1, mu, kr]
+ lineProfile = np.copy(cdf.outint[tStep, 1:cdf.nq[kr]+1, mu, kr])
lineProfile -= lineProfile.min()
lineProfile /= lineProfile.max()
lineProfile *= (y[iwy][0] - y[iwy][-1])
@@ -408,7 +407,7 @@ def make_label(i):
add_legend(ax[1,1])
# Plot 6: Heating, and Core and Wing Contribution Functions
- heat = cdf.bheat1[tStep, iwy]
+ heat = np.copy(cdf.bheat1[tStep, iwy])
if heatPerParticle:
# Only using hydrogen density
heat /= cdf.n1[tStep,:,:6,0].sum(axis=1)[iwy]
@@ -466,19 +465,24 @@ def make_label(i):
hspace=-0.18, wspace=0.01)
if returnData:
- out = {'atomId': cdf.atomid[iel],
+ out = {'atomId': cdf.atomid[0][iel],
'kr': kr,
'iel': iel,
- 'levels': [jTrans, iTrans],
+ 'levels': [jTrans, iTrans],
'labels': [labelI, labelJ],
'emissivity': zTot[np.ix_(iwx, iwy)][:,:, 0, 0],
'opacity': zTot[np.ix_(iwx, iwy)][:,:, 0, 1],
'contFn': zTot[np.ix_(iwx, iwy)][:,:, 1, 1],
+ 'tau': tauq_ny,
'tau1': tau1,
'dVel': dVel,
'xEdges': xEdges,
'yEdges': yEdges,
- 'y': y
+ 'y': y,
+ 'wavelength':wavelength,
+ 'lineProfile':lineProfile,
+ 'iwy':iwy,
+ 'iwx':iwx
}
return out
diff --git a/radynpy/utils/LineProfiles.py b/radynpy/utils/LineProfiles.py
new file mode 100644
index 0000000..128dfc5
--- /dev/null
+++ b/radynpy/utils/LineProfiles.py
@@ -0,0 +1,251 @@
+import numpy as np
+from scipy import constants
+# import sys
+from radynpy.utils.Utils import tradiation as tr
+from radynpy.utils.Utils import tradiation_flux as tr_flux
+
+def profile(cdf, kr, t1 = 0, t2 = 0, trad=False):
+
+ '''
+ A function to return the line intensity, flux, and wavelength
+ of a spectral line transition, from a RADYN simulation
+
+ Parameters
+ __________
+ cdf : the radyn cdf object
+ kr : int
+ the index of the line transition
+ t1 : float
+ the first time at which to extract the line profile
+ default = 0s
+ t2 : float
+ the final time at which to extract the line profile
+ default = final time
+ The output is a time series between t1-->t2
+ trad : boolean
+ set to True to output the line profiles in units of radiation temperature
+ default = false
+
+
+ OUTPUTS: out_dict: A dictionary containing the intensity, flux, and wavelength
+ in units of [erg/s/cm^2/sr/A] (or K), [erg/s/cm^2/A], and [A]
+ Dimensions are:
+ Intensity [time, viewing angle, wavelength]
+ Flux [time, wavelength]
+ Also includes ancillary information such as the
+ nq, the number of wavelength points for the line
+ of interest.
+
+ Orig Written: Graham Kerr, Oct 2020
+
+ Modifications: Graham Kerr, March 2021 :
+ restructing of arrays to have the convention of time
+ being the first dimension.
+ '''
+
+ ########################################################################
+ # Some preliminary set up, constants, and calculation of common terms
+ ########################################################################
+
+ if cdf.cont[kr] != 0:
+ raise ValueError('The transition you entered is not a line transition')
+
+ if t1 !=0:
+ tind1 = np.abs(cdf.time-t1).argmin()
+ t1 = cdf.time[tind1]
+ else:
+ tind1 = 0
+ if t2 !=0:
+ tind2 = np.abs(cdf.time-t1).argmin()
+ else:
+ tind2 = len(cdf.time)-1
+ # t2 = cdf.time[tind2]
+
+
+
+ pi = constants.pi
+ mq = np.nanmax(cdf.nq)
+ # nt = cdf.time.shape[0]
+ nt = len(cdf.time[tind1:tind2+1])
+ ii = cdf.q[0:cdf.nq[kr],kr]
+ wavel = np.zeros(mq,dtype=np.float64)
+ wavel[0:cdf.nq[kr]]=cdf.alamb[kr]/(cdf.q[0:cdf.nq[kr],kr]*cdf.qnorm*1e5/cdf.cc+1)
+ wavel = np.flip(wavel[0:cdf.nq[kr]])
+ line_int = np.zeros([nt,cdf.nmu,mq],dtype=np.float64)
+ line_flux = np.zeros([nt,mq],dtype=np.float64)
+
+ ########################################################################
+ # Extract the line intensity, then sum up the contribution to the line flux
+ ########################################################################
+
+ if trad == True:
+ for ind in range(1,cdf.nq[kr]+1):
+ for tind in range(nt):
+ for muind in range(cdf.nmu):
+ line_int[tind,muind,cdf.nq[kr]-ind] = tr(cdf.outint[tind,ind,muind,kr], wavel[ind-1])
+ line_flux[tind,cdf.nq[kr]-ind] = tr_flux( (2*pi*np.sum(cdf.outint[tind,ind,:,kr]*cdf.wmu*cdf.zmu)), wavel[ind-1])
+ elif trad == False:
+ for ind in range(1,cdf.nq[kr]+1):
+ for tind in range(nt):
+ line_int[tind,:,cdf.nq[kr]-ind] = cdf.outint[tind,ind,:,kr]*cdf.cc*1e8/wavel[ind-1]**2.0
+ line_flux[tind,cdf.nq[kr]-ind] = 2*pi*np.sum(cdf.outint[tind,ind,:,kr]*cdf.wmu*cdf.zmu)*cdf.cc*1e8/wavel[ind-1]**2.0
+
+ rest_wave = cdf.alamb[kr]
+
+ t1 = cdf.time[tind1]
+ t2 = cdf.time[tind2]
+ times = cdf.time[tind1:tind2+1]
+
+ ########################################################################
+ # Output the results
+ ########################################################################
+
+ if trad == False:
+
+ out_dict = {'rest_wave':rest_wave,
+ 'wavelength':wavel,
+ 'line_int':line_int,
+ 'line_flux':line_flux,
+ 't1':t1, 't2':t2,
+ 'tind1':tind1,'tind2':tind2,
+ 'times':times,
+ 'nq':cdf.nq[kr],
+ 'Units':'int in [erg/s/cm^2/sr/A], flux in [erg/s/cm^2/A], wavelength in [A], t1,2 in[s]'}
+
+ elif trad == True:
+
+ out_dict = {'rest_wave':rest_wave,
+ 'wavelength':wavel,
+ 'line_int':line_int,
+ 'line_flux':line_flux,
+ 't1':t1, 't2':t2,
+ 'tind1':tind1,'tind2':tind2,
+ 'times':times,
+ 'nq':cdf.nq[kr],
+ 'Units':'int in [K], flux in [erg/s/cm^2/A], wavelength in [A], t1,2 in[s]'}
+ return out_dict
+
+
+
+def lcurve(cdf, kr, w1 = 0, w2 = 0, t1 = 0, t2 = 0):
+
+ '''
+ A function to return the lightcurves of intensity and flux, for
+ a spectral line transition, from a RADYN simulation.
+
+ Calls LineProfiles.profile to compute the intensity
+ and flux
+
+ Parameters
+ __________
+ cdf : float
+ the radyn readout
+ kr : int
+ the index of the line transition
+ w1 : float
+ the wavelength to start integrating from
+ default = first wavelength in wavelength array
+ w2 : float
+ the wavelength to integrate to
+ default = final wavelength in wavelength array
+ t1 : float
+ the first time at which compute the integrated intensity
+ default = 0s
+ t2 : float
+ the first time at which compute the integrated intensity
+ default = final time
+ The output is a time series from t1 --> t2
+
+ OUTPUTS: out_dict: a dictionary holding the lightcurves of intensity,
+ flux, and the time array, in units of
+ [erg/s/cm^2/sr], [erg/s/cm^2], [s]
+ Dimensions are:
+ Intensity [time, viewing angle]
+ Flux [time]
+ Also includes ancillary information such as the
+ wavelengths and indices over which the integtation
+ was performed.
+
+ Orig Written: Graham Kerr, Oct 2020
+
+ Modifications: Graham Kerr, March 2021 :
+ restructing of arrays to have the convention of time
+ being the first dimension.
+
+ TO DO: Add in option to subtract continuum to return only
+ lightcurves of the line-only intensity or flux
+ '''
+
+ ########################################################################
+ # Some preliminary set up, constants, and calculation of common terms
+ ########################################################################
+
+ if cdf.cont[kr] != 0:
+ raise ValueError('The transition you entered is not a line transition')
+
+ if t1 !=0:
+ tind1 = np.abs(cdf.time-t1).argmin()
+ t1 = cdf.time[tind1]
+ else:
+ tind1 = 0
+ if t2 !=0:
+ tind2 = np.abs(cdf.time-t1).argmin()
+ t2 = cdf.time[tind2]
+ else:
+ tind2 = len(cdf.time)-1
+
+ line = profile(cdf, kr, t1 = t1, t2 = t2)
+
+ if w1 !=0:
+ wind1 = np.abs(line['wavelength']-w1).argmin()
+ w1 = line['wavelength'][wind1]
+ else:
+ wind1 = 0
+ # w1 = line['wavelength'][wind1]
+ if w2 !=0:
+ wind2 = np.abs(line['wavelength']-w2).argmin()
+ w2 = line['wavelength'][wind2]
+ else:
+ wind2 = cdf.nq[kr]-1
+ # w2 = line['wavelength'][wind2]
+
+
+ ########################################################################
+ # Integrate the line intensity and flux over wavelength
+ ########################################################################
+
+
+ lcurve_int = np.zeros([line['line_int'].shape[0],
+ line['line_int'].shape[1]], dtype=np.float64)
+ lcurve_flux = np.zeros([line['line_flux'].shape[0]], dtype=np.float64)
+
+ for tind in range(line['line_int'].shape[0]):
+ lcurve_flux[tind] = np.trapz(line['line_flux'][tind,wind1:wind2+1],line['wavelength'][wind1:wind2+1])
+ for muind in range(line['line_int'].shape[1]):
+ lcurve_int[tind,muind] = np.trapz(line['line_int'][tind,muind,wind1:wind2+1],line['wavelength'][wind1:wind2+1])
+
+ times = cdf.time[tind1:tind2+1]
+
+ ########################################################################
+ # Output the results
+ ########################################################################
+
+ t1 = cdf.time[tind1]
+ t2 = cdf.time[tind2]
+ w1 = line['wavelength'][wind1]
+ w2 = line['wavelength'][wind2]
+
+
+ out_dict = {'lcurve_flux':lcurve_flux,
+ 'lcurve_int':lcurve_int,
+ 'w1':w1, 'w2':w2,
+ 'wind1':wind1,
+ 'wind2':wind2,
+ 't1':t1, 't2':t2,
+ 'tind1':tind1,'tind2':tind2,
+ 'times':times,
+ 'Units':'int in [erg/s/cm^2/sr], flux in [erg/s/cm^2], wavelength in [A], t1,2 in[s]'}
+
+
+
+ return out_dict
diff --git a/radynpy/utils/RadynMovie.py b/radynpy/utils/RadynMovie.py
new file mode 100644
index 0000000..8fac4a7
--- /dev/null
+++ b/radynpy/utils/RadynMovie.py
@@ -0,0 +1,300 @@
+import pandas as pd
+import plotly.express as px
+import plotly.graph_objects as go
+import numpy as np
+
+def rmovie_basicvar(cdf,
+ var = 'tg1',
+ Mm = False,
+ km = False,
+ savefig = False,
+ figname = 'radynvar.html',
+ color = 'steelblue'):
+
+ '''
+ A function to produce an animated figure of RADYN variables.
+ This version is pre-constructed and lets you just input the
+ variable you want to plot. Other variables (such as populations)
+ will require more input, and are separate functions.
+
+ Turns the output into a pandas dataframe, which is then passed to
+ plotly express to create the animated figure
+
+ Parameters
+ __________
+ cdf : The radyn cdf object
+ var : str
+ The variable to plot (default = 'tg1')
+ Mm : Boolean
+ Plot height in Mm (default = False)
+ km : Boolean
+ Plot height in km (default = False)
+ savefig : Boolean
+ Save the figure (html file)
+ figname : str
+ Filename, if saving the output
+
+
+ NOTES :
+ So far, allowed variables are
+ tg1 - temperature
+ ne1 - electron density
+ bheat1 - beam heating rate
+ d1 - mass density
+ vz1 - velocity
+ np - proton density
+
+ Graham Kerr, March 2021
+ '''
+
+ ########################################################################
+ # Some preliminary set up
+ ########################################################################
+
+ if Mm == True:
+ xtitle = 'Height [Mm]'
+ height = cdf.z1/1e8
+ elif km == True:
+ xtitle = 'Height [km]'
+ height = cdf.z1/1e5
+ else:
+ xtitle = 'Height [cm]'
+ height = cdf.z1
+
+ if var == 'tg1':
+ rvar = cdf.tg1
+ ytitle = 'Temperature [K]'
+ ylog = True
+ xlog = False
+ elif var == 'ne1':
+ rvar = cdf.ne1
+ ytitle = 'Electron Density [cm-3]'
+ ylog = True
+ xlog = False
+ elif var == 'bheat1':
+ rvar = cdf.bheat1
+ ytitle = 'Qbeam [erg cm-3 s-1]'
+ ylog = False
+ xlog = False
+ elif var == 'd1':
+ rvar = cdf.d1
+ ytitle = 'Mass Density [g cm-3]'
+ ylog = True
+ xlog = False
+ elif var == 'vz1':
+ rvar = cdf.vz1/1e5
+ ytitle = 'Velocity [km s-1]'
+ ylog = False
+ xlog = False
+ elif var == 'np':
+ rvar = cdf.n1[:,:,5,0]
+ ytitle = 'Proton Density [cm-3]'
+ ylog = True
+ xlog = False
+
+
+
+ template = dict(
+ layout = go.Layout(font = dict(family = "Rockwell", size = 16),
+ title_font = dict(family = "Rockwell", size = 20),
+ plot_bgcolor = 'white',
+ paper_bgcolor = 'white',
+ xaxis = dict(
+ showexponent = 'all',
+ exponentformat = 'e',
+ tickangle = 0,
+ linewidth = 3,
+ showgrid = True,
+ ),
+ yaxis = dict(
+ showexponent = 'all',
+ exponentformat = 'e',
+ linewidth = 3,
+ showgrid = True,
+ anchor = 'free',
+ position = 0,
+ domain = [0.0,1]
+ ),
+ coloraxis_colorbar = dict(
+ thickness = 15,
+ tickformat = '0.2f',
+ ticks = 'outside',
+ titleside = 'right'
+ )
+ ))
+
+ ########################################################################
+ # Build the dataframe
+ ########################################################################
+
+ col1 = ytitle
+ col2 = xtitle
+ time = 'Time [s]'
+ timeind = 'Time index'
+ df_list = []
+ for i in range(len(cdf.time)):
+ data = {col1:rvar[i,:],
+ col2:height[i,:],
+ time: cdf.time[i],
+ timeind: i
+ }
+ df_list.append(pd.DataFrame(data))
+
+ df = pd.concat(df_list)
+
+ ########################################################################
+ # Plot the variable
+ ########################################################################
+
+
+ h1 = 700
+ w1 = 700
+
+ fig1 = px.line(df,
+ x = df.columns[1], y = df.columns[0],
+# animation_group = 'Time [s]',
+ animation_frame = 'Time [s]',
+ log_x = xlog,
+ log_y = ylog,
+ template = template,
+ color_discrete_sequence = [color])
+
+ fig1.show()
+
+ if savefig == True:
+ fig1.write_html(figname)
+
+
+ return df
+
+
+
+
+def rmovie(var1, var2,
+ time = [-10.0],
+ savefig = False,
+ figname = 'radynvar.html',
+ xtitle = 'Var 1',
+ ytitle = 'Var 2',
+ title = ' ',
+ color = 'steelblue',
+ xlog = False, ylog = False):
+
+ '''
+ A function to produce an animated figure of RADYN variables.
+ This version is 'dumb' and just plots col1 vs col2 without any
+ axes labels, unless passed through the function fall.
+
+ Variables must be input as [time, dim1]
+
+ Turns the output into a pandas dataframe, which is then passed to
+ plotly express to create the animated figure
+
+ Parameters
+ __________
+ var1 : float
+ The variable to plot on the x-axis [time, dim1]
+ var2 : float
+ The variable to plot on the y-axis [time, dim1]
+ xtitle : str
+ The xaxis label (default "Var 1")
+ ytitle : str
+ The xaxis label (default "Var 2")
+ title : str
+ A plot title (default " ")
+ savefig : Boolean
+ Save the figure (html file)
+ figname : str
+ Filename, if saving the output
+ xlog : boolean
+ Default is false. Set to True to have log x-axis
+ ylog : boolean
+ Default is false. Set to True to have log y-axis
+
+ NOTES :
+
+
+ Graham Kerr, March 2021
+ '''
+
+ ########################################################################
+ # Some preliminary set up
+ ########################################################################
+
+ if time[0] == -10:
+ time = np.arange(0,var1.shape[0])
+ col3 = 'Time [index]'
+ else:
+ col3 = 'Time [s]'
+
+ template = dict(
+ layout = go.Layout(font = dict(family = "Rockwell", size = 16),
+ title_font = dict(family = "Rockwell", size = 20),
+ plot_bgcolor = 'white',
+ paper_bgcolor = 'white',
+ xaxis = dict(
+ showexponent = 'all',
+ exponentformat = 'e',
+ tickangle = 0,
+ linewidth = 3,
+ showgrid = True,
+ ),
+ yaxis = dict(
+ showexponent = 'all',
+ exponentformat = 'e',
+ linewidth = 3,
+ showgrid = True,
+ anchor = 'free',
+ position = 0,
+ domain = [0.0,1]
+ ),
+ coloraxis_colorbar = dict(
+ thickness = 15,
+ tickformat = '0.2f',
+ ticks = 'outside',
+ titleside = 'right'
+ )
+ ))
+
+ ########################################################################
+ # Build the dataframe
+ ########################################################################
+
+ col1 = xtitle
+ col2 = ytitle
+ df_list = []
+ for i in range(len(time)):
+ data = {col1:var1[i,:],
+ col2:var2[i,:],
+ col3: time[i],
+ }
+ df_list.append(pd.DataFrame(data))
+
+ df = pd.concat(df_list)
+
+
+ ########################################################################
+ # Plot the variable
+ ########################################################################
+
+
+ h1 = 700
+ w1 = 700
+
+ fig1 = px.line(df,
+ x = df.columns[0], y = df.columns[1],
+# animation_group = 'Time [s]',
+ animation_frame = df.columns[2],
+ log_x = xlog,
+ log_y = ylog,
+ title = title,
+ color_discrete_sequence = [color],
+ template = template)
+
+ fig1.show()
+
+ if savefig == True:
+ fig1.write_html(figname)
+
+
+ return df
\ No newline at end of file
diff --git a/radynpy/utils/Utils.py b/radynpy/utils/Utils.py
index b9ab6bb..cb6976b 100644
--- a/radynpy/utils/Utils.py
+++ b/radynpy/utils/Utils.py
@@ -264,3 +264,508 @@ def planck_fn(wvls, tg):
bbflux *= 1e-8 / np.pi
return bbflux.squeeze()
+
+
+def prfhbf_rad(wvls = [], Z = 1, n=6, *args):
+
+ '''
+ A function to return the absorption
+ cross section fr a hydrogenic ion
+ *** Originaly from abshyd.pro subroutine
+ in RADYN, and uses Fortran indexing.
+ Using this code to ensure that we
+ compute the upper level contribution
+ to the opacity in the same way that
+ RADYN does internally.
+ Parameters
+ __________
+ wvls : float
+ the wavelength(s)
+ n : int, optional
+ The level from which absorption
+ takes place. FORTRAN NUMBERING, legacy
+ from RADYN (default = 6, i.e. level 6)
+ Z : int, optional
+ charge (default = 1)
+ Graham Kerr, Feb 19th 2020
+ '''
+ wvls = np.array(wvls)
+
+ prfhbf = np.zeros([len(wvls)],dtype=float)
+
+ wl0 = 911.7535278/(Z*Z)*n*n
+
+ for i in range(len(wvls)):
+
+ if (wvls[i] > wl0):
+
+ prfhbf[i] = 0.0
+
+ else:
+
+ frq = 2.9979e18/wvls[i]
+
+ gau = gaunt_factor(n,frq)
+
+ pr0 = 1.04476e-14*Z*Z*Z*Z
+
+ a5 = n**5
+ a5 = np.float(a5)
+
+ wm = wvls[i]*1.0e-4
+
+ wm3 = wm*wm*wm
+
+ prfhbf[i] = pr0*wm3*gau/a5
+
+
+ return prfhbf
+
+def hminus_pops(tg1, ne1, nh, bhmin=1.0, *args):
+
+ '''
+ A function to return the H minus populations
+ Parameters
+ __________
+ tg1 : float
+ the temperature stratification
+ nh : float
+ the population densities of hydrogen
+ bhmin : float, optional
+ the departure coefficient from LTE
+ for H minus (default = 0).
+ Make sure the tg, and nh have dimensions of
+ [timesteps, height], and [timesteps, height, level]
+
+ The density of H minus from
+ Vernazza et al. 1976, ApJS 30
+ Graham Kerr, Feb 19th 2020
+ '''
+ ndep = (tg1.shape)[1]
+ num_times = (tg1.shape)[0]
+ num_levs = (nh.shape)[2]
+
+ nhmin = np.zeros([num_times,ndep],dtype=float)
+ totnhi = np.sum(nh[:,:,0:num_levs-1],axis = 2) # Density of neutral hydrogen
+
+ for k in range(num_times):
+ nhmin[k,:] = (
+ 1.0354e-16 * bhmin * ne1[k,:] * totnhi[k,:] *
+ tg1[k,:]**(-1.5) * np.exp(8762e0/tg1[k,:])
+ )
+
+ return nhmin
+
+def transp_scat(tau=[],x=[],sc=[],scat=[],
+ *args):
+
+ '''
+ A function to return the average intensity J at some
+ wavelength
+ Parameters
+ __________
+
+ tau : float
+ the optical depth at some standard wavelength
+ (usually 5000A)
+ x : float
+ the ratio of total opacity to the opacity at some
+ standard wavelength
+ sc : float
+ epsilon B
+ scat: float
+ 1 - epsilon
+
+ Graham Kerr, Feb 20th 2020
+ '''
+ ntimes = (tau.shape)[0]
+ ndep = (tau.shape)[1]
+ nwave = (x.shape)[0]
+
+
+ nmu=3
+ wmu=np.array([0.277778,0.444444,0.277778], dtype=float)
+ xmu=np.array([0.112702,0.500000,0.887298], dtype=float)
+
+ itran=0
+
+ sp1=np.zeros([nmu,nmu,nwave,ndep,ntimes],dtype=float)
+ sp2=np.zeros([nmu,nmu,nwave,ndep,ntimes],dtype=float)
+ sp3=np.zeros([nmu,nmu,nwave,ndep,ntimes],dtype=float)
+ a1=np.zeros([nwave,ndep,ntimes],dtype=float)
+ c1=np.zeros([nwave,ndep,ntimes],dtype=float)
+ p=np.zeros([nmu,nwave,ndep,ntimes],dtype=float)
+ pms=np.zeros([nwave,ndep,ntimes],dtype=float)
+ tauq=np.zeros([nwave,ndep,ntimes],dtype=float)
+ dtauq=np.zeros([nwave,ndep,ntimes],dtype=float)
+ jnu=np.zeros([nwave,ndep,ntimes],dtype=float)
+
+
+ cmu=0.5/xmu[0]
+ for k in range(ntimes):
+ for i in range(nwave):
+ #
+ # j=1: upper boundary
+ #
+ dtauq[i,1,k]=((x[i,0,k]+x[i,1,k])*(tau[k,1]-tau[k,0]))*cmu
+ a1[i,0,k]=1./dtauq[i,1,k]
+ t = tau[k,0]*x[i,0,k]*2.0*cmu
+ tauq[i,0,k]=t
+ dtauq[i,0,k]=t
+ dtauq[i,1:ndep,k]=(x[i,1:ndep,k]+x[i,0:ndep-1,k])*(tau[k,1:ndep]-tau[k,0:ndep-1])*cmu
+ for j in range(1,ndep):
+ tauq[i,j,k]=tauq[i,j-1,k]+dtauq[i,j,k]
+ #
+ # calculate a1 and c1
+ #
+ a1[i,1:ndep-1,k] = 2./(dtauq[i,1:ndep-1,k]+dtauq[i,2:ndep,k])/dtauq[i,1:ndep-1,k]
+ c1[i,1:ndep-1,k] = 2./(dtauq[i,1:ndep-1,k]+dtauq[i,2:ndep,k])/dtauq[i,2:ndep,k]
+
+ for mu in range(nmu):
+ xmu1=xmu[mu]/xmu[0]
+ xmu2=xmu1*xmu1
+
+ #
+ # k=1: upper boundary
+ #
+ t=tauq[i,0,k]/xmu1
+ if (t < 0.01):
+ ex1=t*(1.-t*(0.5-t*(0.1666667-t*0.041666667)))
+ elif(t < 20.):
+ ex1=1.-np.exp(-t)
+ else:
+ ex1=1.
+
+ ex = 1.-ex1
+ sp1[mu,mu,i,0,k] = 0.
+ sp2[mu,mu,i,0,k] = 1.0+2.0*xmu1*a1[i,0,k]
+ sp3[mu,mu,i,0,k] = -2.0*xmu2*a1[i,0,k]*a1[i,0,k]
+ fact = 1.0+2.0*xmu1*a1[i,0,k]*ex1
+ sp2[mu,mu,i,0,k] = sp2[mu,mu,i,0,k]/fact
+ sp3[mu,mu,i,0,k] = sp3[mu,mu,i,0,k]/fact
+ p[mu,i,0,k]=sc[i,0,k]
+
+ #
+ # interior
+ #
+ sp1[mu,mu,i,1:ndep-1,k]=-xmu2*a1[i,1:ndep-1,k]
+ sp2[mu,mu,i,1:ndep-1,k]=1.0
+ sp3[mu,mu,i,1:ndep-1,k]=-xmu2*c1[i,1:ndep-1,k]
+ p[mu,i,1:ndep-1,k]=sc[i,1:ndep-1,k]
+
+ #
+ # k=ndep: lower boundary
+ #
+ sp1[mu,mu,i,ndep-1,k]= -1.0
+ sp2[mu,mu,i,ndep-1,k]=dtauq[i,ndep-1,k]/xmu1+0.5*dtauq[i,ndep-1,k]**2/xmu2
+ sp3[mu,mu,i,ndep-1,k]=0.0
+ sk=(dtauq[i,ndep-1,k]/xmu1+0.5*dtauq[i,ndep-1,k]**2/xmu2)+1.0
+ p[mu,i,ndep-1,k]=sc[i,ndep-1,k]*sk - sc[i,ndep-2,k]
+
+ #
+ # non-diagonal elements
+ #
+ for mu2 in range(nmu):
+ sp2[mu,mu2,i,0,k] = sp2[mu,mu2,i,0,k]-scat[i,0,k]*wmu[mu2]
+ sp2[mu,mu2,i,1:ndep-1,k] = sp2[mu,mu2,i,1:ndep-1,k]-scat[i,1:ndep-1,k]*wmu[mu2]
+ sp2[mu,mu2,i,ndep-1,k] = ( sp2[mu,mu2,i,ndep-1,k]-scat[i,ndep-1,k]*wmu[mu2]*sk +
+ scat[i,ndep-2,k]*wmu[mu2])
+ sp1[mu,mu2,i,ndep-1,k]= sp1[mu,mu2,i,ndep-1,k]+scat[i,ndep-2,k]*wmu[mu2]
+
+ #
+ # eliminate subdiagonal
+ #
+ for j in range(0,ndep-1):
+ sp1p = sp1[:,:,i,j+1,k]
+ sp2k = sp2[:,:,i,j,k]
+ sp3k = sp3[:,:,i,j,k]
+ f = -np.matmul(sp1p,np.linalg.inv(sp2k-sp3k))
+ #f = -sp1p
+ p[:,i,j+1,k] = p[:,i,j+1,k]+np.matmul(f,(p[:,i,j,k]))
+ sp2[:,:,i,j+1,k] = sp2[:,:,i,j+1,k]+np.matmul(f,sp2k)
+ sp2[:,:,i,j,k] = sp2[:,:,i,j,k] - sp3[:,:,i,j,k]
+
+ sp2[:,:,i,ndep-1,k]=sp2[:,:,i,ndep-1,k]-sp3[:,:,i,ndep-1,k]
+
+ #
+ # backsubstitute
+ #
+ p[:,i,ndep-1,k] = np.matmul( np.linalg.inv(sp2[:,:,i,ndep-1,k]) ,p[:,i,ndep-1,k])
+
+
+ for j in range(ndep-2, -1, -1):
+ tmp = np.linalg.inv(sp2[:,:,i,j,k])
+ tmp2 = np.matmul(sp3[:,:,i,j,k],p[:,i,j+1,k])
+ tmp3 = p[:,i,j,k] - tmp2
+ p[:,i,j,k]= np.matmul( tmp, tmp3 )
+
+
+ for mu in range(nmu):
+ jnu[i,:,k] = jnu[i,:,k] + wmu[mu]*p[mu,i,:,k]
+
+
+ return jnu
+
+
+def wnstark2(n=[],ne=[],n1=[],tg=[],
+ *args):
+
+ '''
+ The occupational probability.
+ That is, the probability of an atom being perturbed by a net electric
+ microfield that is below the critical field strength that would
+ dissolve/destroy level n. Used for Landau-Zener effects
+ (see Kowalski et al 2015)
+ Parameters
+ __________
+
+ n : float
+ pseudo quantum number
+ ne : float
+ the electron density in height
+ n1 : float
+ the ground state H density in height
+ tg: float
+ the gas temperature in height
+ Graham Kerr, Feb 27th 2020
+ '''
+
+ Zjk = 1.0e0
+ Zr = 1.0e0
+ el = 4.803e-10
+
+ a0 = 5.2981e-9
+ F0 = 1.25e-9 * Zjk * ne**(2./3.)
+
+ betacrit = (2.0e0 * n + 1.0e0) / (6.0e0 * n**4.0 * (n+1.0e0)**2) * el / (a0**2.0 * F0)
+
+ if n < 3:
+ Kn = 1.0e0
+ else:
+ Kn = 16.0e0/3.0e0 * n / (n+1.0e0)**2.0
+
+ P1 = 0.1402e0
+ P2 = 0.1285e0
+ P3 = 1e0
+ P4 = 3.15e0
+ P5 = 4.0e0
+
+ a = 0.09e0 * (ne**(1.0e0/6.0e0)) / tg**0.5e0
+
+ X = (1.0e0 + P3*a)**P4
+ C1 = P1*(X + P5*Zr*a**3)
+ C2 = P2*X
+
+ lilf = C1 * betacrit**3 / (1.0e0 + C2 * betacrit**(1.5e0))
+
+ Q = lilf / (1.0e0+lilf)
+
+ wn_neutral = np.exp(-1.0e0 * 4.0e0 * np.pi/3.0e0 * n1 * (n**2.0e0*a0+a0)**3.0e0)
+
+ return wn_neutral*Q
+
+
+
+def poph2(cdf,*args):
+
+ '''
+ The population of H2 molecules from LTE dissociation
+ equilibrium
+ Parameters
+ __________
+
+ cdf :
+ radynpy cdf output
+
+ Graham Kerr, Feb 27th 2020
+ **** Could probably re-write to only input certain
+ variables but for now they are all passed via
+ the cdf file
+ '''
+ grph = 2.2690989266985202e-24
+ nt = (cdf.time.shape)[0]
+ nh21t = np.zeros([nt, cdf.ndep])
+
+ for k in range(nt):
+ for j in range(cdf.ndep):
+
+ xkh2k,dxkh2k,eh2k,deh2k = molfys(cdf.tg1[k,j])
+ a1 = 8.e0*cdf.d1[k,j]/grph/xkh2k
+
+ if (a1 > 1.e-03):
+ nh21t[k,j]=(4.e0*cdf.d1[k,j]/grph+xkh2k-xkh2k*np.sqrt(1.e0+a1))/8.e0
+ else:
+ a2=(cdf.d1[k,j]/grph)**2/xkh2k
+ a3=(cdf.d1[k,j]/grph)**3*4.e0/xkh2k/xkh2k
+ nh21t[k,j]=a2-a3
+
+ return nh21t
+
+
+
+
+
+def molfys(tg,*args):
+
+ '''
+ Gives the dissociation constant xkh2 (=n(h i)*n(h i)/n(h2) )
+ expressed in number per cm3, and the sum of dissociation,
+ rotation and vibration energies, deh2 for h2 (expressed in
+ ergs per molecule).
+ The data are from vardya, m.n.r.a.s. 129, 205 (1965) and earlier
+ references. the dissociation constant for h2 is from tsuji,
+ astron. astrophys. 1973.
+ The logarithmic temperature derivatives are also returned in
+ dxkh2 and deh2
+ Parameters
+ __________
+
+ tg : float
+ the gas temperature
+
+ Graham Kerr, Feb 27th 2020
+ '''
+ a1 = np.array([12.739e0,-5.1172e0,1.2572e-1,-1.4149e-2,6.3021e-4])
+ b1 = np.array([2.6757e0,-1.4772e0,0.60602e0,-0.12427e0,0.009750e0])
+ bk=1.380662e-16
+ ee=1.6021890e-12
+
+ te = np.zeros([5],dtype=float)
+ tes = np.zeros([7],dtype=float)
+ tex=5040.e0/tg
+ te[0]=1.e0
+ for k in range(1,5):
+ te[k]=te[k-1]*tex
+ tes[2:7] = te
+ tes[1] = 1.e0/tex
+ tes[0] = te[1]/tex
+
+ # dissociation constant for h2
+ xkh2=0.e0
+ for k in range(5):
+ xkh2=a1[k]*te[k]+xkh2
+ xkh2=10.e0**xkh2
+ xkh2=xkh2/bk/tg
+
+ # logarithmic temperature derivative of dissociation constant for h2
+ dxkh2=0.e0
+ for k in range(1,5):
+ dxkh2=k*a1[k]*te[k-1]+dxkh2
+ dxkh2=-dxkh2*np.log(10.e0)*xkh2*tex
+ dxkh2=dxkh2-xkh2
+
+ # internal energy per molecule
+ eh2=0.e0
+ for k in range(5):
+ eh2=b1[k]*tes[2+k-1]+eh2
+ eh2=eh2*8.617e-5*5040.e0-4.476e0
+ eh2=eh2*ee
+
+ # logarithmic temperature derivative of internal energy
+
+ deh2=0.e0
+ for k in range(5):
+ deh2=np.float(k-1)*b1[k]*tes[2+k-2]+deh2
+ deh2=-deh2*8.617e-5*5040.e0*tex
+ deh2=deh2*ee
+
+ return xkh2,dxkh2,eh2,deh2
+
+
+
+def closest_ind(x, y=[], *args):
+ '''
+ Finds the index in an array X such that X(ind) is
+ closest to y. Returns the index and the difference
+ Parameters
+ __________
+
+ x :
+ an array
+ y :
+ the value(s) within x to find
+
+ Graham Kerr, March 2nd 2020
+ '''
+ numvals = len(y)
+
+ inds = np.zeros([numvals],dtype=int)
+
+ for i in range(numvals):
+ inds[i] = (np.abs(x - y[i])).argmin()
+
+ diffs = x[inds] - y
+
+ return inds, diffs
+
+
+
+def tradiation(intens, lamb):
+ '''
+ Graham Kerr
+ NASA/GSFC
+ 9th June 2021
+
+ NAME: tradiation.py
+
+ PURPOSE: Convert from cgs intensity to the radiation temperature.
+ That is, what would be the temperature of a blackbody that
+ emits at the given intensity?
+
+ INPUTS: intens -- float array
+ the intensity in units of erg/s/cm^2/sr/Hz
+ lamb -- float array
+ the wavelengh array in units of angstroms
+
+ OPTIONAL
+ INPUTS:
+
+
+ OUTPUTS: trad -- the intensity in units of K
+
+
+ NOTES:
+ '''
+
+ c = 2.99792458e10
+ h = 6.626176e-27
+ k = 1.380662e-16
+ l = lamb*1e-8
+ tRad = h * c / k / l / np.log(2.0 * h * c / intens / (l**3) + 1.0)
+
+ return tRad
+
+def tradiation_flux(flux, lamb):
+ '''
+ Graham Kerr
+ NASA/GSFC
+ 9th June 2021
+
+ NAME: tradiation.py
+
+ PURPOSE: Convert from cgs flux to the radiation temperature.
+ That is, what would be the temperature of a blackbody that
+ emits at the given intensity?
+
+ INPUTS: flux -- float array
+ the flux in units of erg/s/cm^2/Hz
+ lamb -- float array
+ the wavelengh array in units of angstroms
+
+ OPTIONAL
+ INPUTS:
+
+
+ OUTPUTS: trad -- the intensity in units of K
+
+
+ NOTES:
+ '''
+
+ c = 2.99792458e10
+ h = 6.626176e-27
+ k = 1.380662e-16
+ l = lamb*1e-8
+ tRad_flux = h * c / k / l / np.log(2.0 * np.pi * h * c / flux / (l**3) + 1.0)
+
+ return tRad_flux
diff --git a/setup.py b/setup.py
index 2608240..c09af1b 100644
--- a/setup.py
+++ b/setup.py
@@ -13,7 +13,7 @@ def readme():
author_email='c.osborne.1@research.gla.ac.uk',
license='MIT',
packages=find_packages(),
- install_requires=['numpy', 'scikit-image', 'matplotlib', 'scipy', 'colour', 'palettable', 'cdflib'],
+ install_requires=['numpy', 'scikit-image', 'matplotlib', 'scipy', 'colour', 'palettable', 'cdflib', 'plotly','pandas'],
classifiers=[
'Development Status :: 3 - Alpha',
'Intended Audience :: Science/Research',