diff --git a/tests/manual_check/comparison/comparison_getE.py b/tests/manual_check/comparison/comparison_getE.py deleted file mode 100644 index ff358e5c..00000000 --- a/tests/manual_check/comparison/comparison_getE.py +++ /dev/null @@ -1,28 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -from exojax.dynamics.getE import getE -# make example by PyAstronomy -from PyAstronomy.pyasl.asl import MarkleyKESolver as MKE - -m = MKE() -print(m.getE(M=0.5, e=0.3)) -marr = np.linspace(0.0, 4*np.pi, 1000) -ea = [] -for meach in marr: - ea.append(m.getE(M=meach, e=0.3)) -ea = np.array(ea) - -fig = plt.figure() -ax = fig.add_subplot(211) -plt.plot(marr, ea) -plt.plot(marr, getE(marr, 0.3)) - -ax.set_ylabel('Eccentric anomary') - -ax = fig.add_subplot(212) -plt.plot(marr, ea-getE(marr, 0.3)) - -plt.xlabel('Mean anomary') -ax.set_ylabel('Difference') -plt.savefig('getE.png') -plt.show() diff --git a/tests/manual_check/comparison/dit_hitran_CO.py b/tests/manual_check/comparison/dit_hitran_CO.py deleted file mode 100644 index 08895bbe..00000000 --- a/tests/manual_check/comparison/dit_hitran_CO.py +++ /dev/null @@ -1,97 +0,0 @@ -from jax import config -from exojax.spec.hitran import normalized_doppler_sigma -from exojax.spec.set_ditgrid import ditgrid_log_interval -from exojax.spec import api -from exojax.spec.hitran import line_strength, doppler_sigma, gamma_hitran, gamma_natural -from exojax.spec.lpf import auto_xsection as lpf_xsection -from exojax.spec import initspec -from exojax.spec.dit import xsvector as dit_xsvector -from exojax.spec.lpf import xsvector as lpf_xsvector -from exojax.spec import make_numatrix0 -import numpy as np -import tqdm -import jax.numpy as jnp -from jax import vmap -import seaborn as sns -import matplotlib.pyplot as plt -plt.style.use('bmh') - -nus = np.linspace((4000), (4500.0), 3000000, dtype=np.float64) -mdbCO = api.MdbHitran('~/exojax/data/CO/05_hit12.par', nus) - -Mmol = 28.010446441149536 -Tref = 296.0 -Tfix = 1000.0 -Pfix = 1.e-3 - -# USE TIPS partition function -Q296 = np.array([107.25937215917970, 224.38496958496091, 112.61710362499998, - 660.22969049609367, 236.14433662109374, 1382.8672147421873]) -Q1000 = np.array([382.19096582031250, 802.30952197265628, 402.80326733398437, - 2357.1041210937501, 847.84866308593757, 4928.7215078125000]) -qr = Q1000/Q296 - -qt = np.ones_like(mdbCO.isoid, dtype=np.float64) -for idx, iso in enumerate(mdbCO.uniqiso): - mask = mdbCO.isoid == iso - qt[mask] = qr[idx] - -Sij = line_strength(Tfix, mdbCO.logsij0, mdbCO.nu_lines, mdbCO.elower, qt, mdbCO.Tref) -gammaL = gamma_hitran(Pfix, Tfix, Pfix, mdbCO.n_air, - mdbCO.gamma_air, mdbCO.gamma_self) -# + gamma_natural(A) #uncomment if you inclide a natural width -sigmaD = doppler_sigma(mdbCO.nu_lines, Tfix, Mmol) - -cnu, indexnu, pmarray = initspec.init_dit(mdbCO.nu_lines, nus) -sigmaD_grid = ditgrid_log_interval(sigmaD, dit_grid_resolution=0.1) -gammaL_grid = ditgrid_log_interval(gammaL, dit_grid_resolution=0.1) - -xs_dit_lp = dit_xsvector(cnu, indexnu, pmarray, sigmaD, - gammaL, Sij, nus, sigmaD_grid, gammaL_grid) -wls_dit = 100000000/nus - -#ref (direct) -d = 10000 -ll = mdbCO.nu_lines -xsv_lpf_lp = lpf_xsection(nus, ll, sigmaD, gammaL, Sij, memory_size=30) - -config.update('jax_enable_x64', True) -xs_dit_lp_f64 = dit_xsvector( - cnu, indexnu, pmarray, sigmaD, gammaL, Sij, nus, sigmaD_grid, gammaL_grid) - - -# PLOT -llow = 2300.4 -lhigh = 2300.7 -tip = 20.0 -fig = plt.figure(figsize=(12, 3)) -ax = plt.subplot2grid((12, 1), (0, 0), rowspan=8) -plt.plot(wls_dit, xsv_lpf_lp, label='Direct', - color='C0', markersize=3, alpha=0.99) -plt.plot(wls_dit, xs_dit_lp, ls='dashed', color='C1', alpha=0.7, label='DIT') - -plt.ylim(1.1e-28, 1.e-17) -# plt.ylim(1.e-27,3.e-20) -plt.yscale('log') - -plt.xlim(llow*10.0-tip, lhigh*10.0+tip) -plt.legend(loc='upper right') -plt.ylabel(' cross section', fontsize=10) -#plt.text(22986,3.e-21,"$P=10^{-3}$ bar") -plt.xlabel('wavelength [$\AA$]') - -ax = plt.subplot2grid((12, 1), (8, 0), rowspan=4) -plt.plot(wls_dit, np.abs(xs_dit_lp/xsv_lpf_lp-1.)*100, - lw=1, alpha=0.5, color='C1', label='MODIT (F32)') -plt.plot(wls_dit, np.abs(xs_dit_lp_f64/xsv_lpf_lp-1.)*100, - lw=1, alpha=1, color='C2', label='MODIT (F64)') -plt.yscale('log') -plt.ylabel('difference (%)', fontsize=10) -plt.xlim(llow*10.0-tip, lhigh*10.0+tip) -plt.ylim(0.01, 100.0) -plt.xlabel('wavelength [$\AA$]') -plt.legend(loc='upper left') - -plt.savefig('comparison_dit.png', bbox_inches='tight', pad_inches=0.0) -plt.savefig('comparison_dit.pdf', bbox_inches='tight', pad_inches=0.0) -plt.show() diff --git a/tests/manual_check/comparison/hi_dynamic_range/errCO_DIT.py b/tests/manual_check/comparison/hi_dynamic_range/errCO_DIT.py deleted file mode 100644 index 312bf9eb..00000000 --- a/tests/manual_check/comparison/hi_dynamic_range/errCO_DIT.py +++ /dev/null @@ -1,118 +0,0 @@ -import numpy as np -import tqdm -import jax.numpy as jnp -from jax import vmap -import seaborn as sns -import matplotlib.pyplot as plt -plt.style.use('bmh') -from exojax.spec import make_numatrix0 -from exojax.spec.lpf import xsvector as lpf_xsvector -from exojax.spec.dit import xsvector as dit_xsvector -from exojax.spec import initspec -from exojax.spec.lpf import auto_xsection as lpf_xsection -from exojax.spec.hitran import line_strength, doppler_sigma, gamma_hitran, gamma_natural -from exojax.spec import api -from exojax.spec.set_ditgrid import ditgrid_log_interval -from exojax.spec.hitran import normalized_doppler_sigma - - -def comperr(Nnu,plotfig=False): - - nus=np.linspace(1.e7/2700,1.e7/2100.,Nnu,dtype=np.float64) - - mdbCO=api.MdbHitran('~/exojax/data/CO/05_hit12.par',nus) - - Mmol=28.010446441149536 - Tref=296.0 - Tfix=1000.0 - Pfix=1.e-3 # - - #USE TIPS partition function - Q296=np.array([107.25937215917970,224.38496958496091,112.61710362499998,\ - 660.22969049609367,236.14433662109374,1382.8672147421873]) - Q1000=np.array([382.19096582031250,802.30952197265628,402.80326733398437,\ - 2357.1041210937501,847.84866308593757,4928.7215078125000]) - qr=Q1000/Q296 - - qt=np.ones_like(mdbCO.isoid,dtype=np.float64) - for idx,iso in enumerate(mdbCO.uniqiso): - mask=mdbCO.isoid==iso - qt[mask]=qr[idx] - - Sij=line_strength(Tfix,mdbCO.logsij0,mdbCO.nu_lines,mdbCO.elower,qt,mdbCO.Tref) - gammaL = gamma_hitran(Pfix,Tfix,Pfix, mdbCO.n_air, mdbCO.gamma_air, mdbCO.gamma_self) - #+ gamma_natural(A) #uncomment if you inclide a natural width - sigmaD=doppler_sigma(mdbCO.nu_lines,Tfix,Mmol) - sigmaD_grid=ditgrid_log_interval(sigmaD) - gammaL_grid=ditgrid_log_interval(gammaL) - - - cnu,indexnu,dLarray=initspec.init_dit(mdbCO.nu_lines,nus) - - xs_dit_lp=dit_xsvector(cnu,indexnu,dLarray,sigmaD,gammaL,Sij,nus,sigmaD_grid,gammaL_grid) - wls_dit = 100000000/nus - - #ref (direct) - d=10000 - ll=mdbCO.nu_lines - xsv_lpf_lp=lpf_xsection(nus,ll,sigmaD,gammaL,Sij,memory_size=30) - - - dif=xs_dit_lp/xsv_lpf_lp-1. - med=np.median(dif) - iju=22940. - ijd=26400. - limu,limd=np.searchsorted(wls_dit[::-1],[iju,ijd]) - std=np.std(dif[::-1][limu:limd]) - - return med,std,ijd,iju,wls_dit,xs_dit_lp,xsv_lpf_lp,dif - - -if __name__=="__main__": - import matplotlib - m,std,ijd,iju,wls_dit,xs_dit_lp,xsv_lpf_lp,dif=comperr(200000) - m1,std1,ijd1,iju1,wls_dit1,xs_dit_lp1,xsv_lpf_lp1,dif1=comperr(3000000) - - #PLOT - plotfig=True - if plotfig: - matplotlib.rcParams['agg.path.chunksize'] = 100000 - llow=2300.4 - lhigh=2300.7 - tip=2.0 - fig=plt.figure(figsize=(12,3)) - ax=plt.subplot2grid((12, 1), (0, 0),rowspan=8) - plt.plot(wls_dit1,xsv_lpf_lp1,label="Direct",color="C0",alpha=0.3,markersize=3) - plt.plot(wls_dit,xs_dit_lp,color="C1",lw=1,alpha=0.3) - plt.plot(wls_dit1,xs_dit_lp1,color="C2",lw=1,alpha=0.5,ls="dashed") - - plt.xlim(ijd,iju) -# plt.xlim(llow*10-tip,lhigh*10+tip) - - plt.ylim(1.1e-35,1.e-17) -# plt.ylim(1.e-27,3.e-17) - plt.yscale("log") - -# plt.xlim(llow*10.0-tip,lhigh*10.0+tip) - plt.legend(loc="upper right") - plt.ylabel(" cross section",fontsize=10) - #plt.text(22986,3.e-21,"$P=10^{-3}$ bar") - plt.xlabel('wavelength [$\AA$]') - - ax=plt.subplot2grid((12, 1), (8, 0),rowspan=4) - plt.plot(wls_dit,np.abs((dif)*100),alpha=0.2,color="C1") - plt.plot(wls_dit1,np.abs((dif1)*100),alpha=0.5,color="C2") - - plt.ylabel("difference (%)",fontsize=10) - plt.xlim(ijd,iju) -# plt.xlim(llow*10-tip,lhigh*10+tip) - - plt.ylim(0.1,10000.0) -# plt.ylim(-10*100*std,10*100*std) - plt.yscale("log") - plt.xlabel('wavelength [$\AA$]') - plt.legend(loc="upper left") - - plt.savefig("fig/comparison_dit.png", bbox_inches="tight", pad_inches=0.0) - plt.savefig("fig/comparison_dit.pdf", bbox_inches="tight", pad_inches=0.0) - plt.show() diff --git a/tests/manual_check/comparison/hi_dynamic_range/errCO_MODIT.py b/tests/manual_check/comparison/hi_dynamic_range/errCO_MODIT.py deleted file mode 100644 index 4237c1fd..00000000 --- a/tests/manual_check/comparison/hi_dynamic_range/errCO_MODIT.py +++ /dev/null @@ -1,128 +0,0 @@ -import numpy as np -import tqdm -import jax.numpy as jnp -from jax import vmap -import seaborn as sns -import matplotlib.pyplot as plt -plt.style.use('bmh') -from exojax.spec import make_numatrix0 -from exojax.spec.lpf import xsvector as lpf_xsvector -from exojax.spec.modit import xsvector as modit_xsvector -from exojax.spec import initspec -from exojax.spec.lpf import auto_xsection as lpf_xsection -from exojax.spec.hitran import line_strength, doppler_sigma, gamma_hitran, gamma_natural -from exojax.spec import api -from exojax.spec.set_ditgrid import ditgrid_log_interval -from exojax.spec.hitran import normalized_doppler_sigma - -#F64/F32 -from jax import config -config.update("jax_enable_x64", True) - -def comperr(Nnu,plotfig=False): - - nus=np.logspace(np.log10(1.e7/2700),np.log10(1.e7/2100.),Nnu,dtype=np.float64) - -# nus=np.logspace(np.log10(3000),np.log10(6000.0),Nnu,dtype=np.float64) - mdbCO=api.MdbHitran('~/exojax/data/CO/05_hit12.par',nus) - - Mmol=28.010446441149536 - Tref=296.0 - Tfix=1000.0 - Pfix=1.e-3 # - - #USE TIPS partition function - Q296=np.array([107.25937215917970,224.38496958496091,112.61710362499998,\ - 660.22969049609367,236.14433662109374,1382.8672147421873]) - Q1000=np.array([382.19096582031250,802.30952197265628,402.80326733398437,\ - 2357.1041210937501,847.84866308593757,4928.7215078125000]) - qr=Q1000/Q296 - - qt=np.ones_like(mdbCO.isoid,dtype=np.float64) - for idx,iso in enumerate(mdbCO.uniqiso): - mask=mdbCO.isoid==iso - qt[mask]=qr[idx] - - Sij=line_strength(Tfix,mdbCO.logsij0,mdbCO.nu_lines,mdbCO.elower,qt,mdbCO.Tref) - gammaL = gamma_hitran(Pfix,Tfix,Pfix, mdbCO.n_air, mdbCO.gamma_air, mdbCO.gamma_self) - #+ gamma_natural(A) #uncomment if you inclide a natural width - sigmaD=doppler_sigma(mdbCO.nu_lines,Tfix,Mmol) - - cnu,indexnu,R,dq=initspec.init_modit(mdbCO.nu_lines,nus) - nsigmaD=normalized_doppler_sigma(Tfix,Mmol,R) - ngammaL=gammaL/(mdbCO.nu_lines/R) - ngammaL_grid = ditgrid_log_interval(ngammaL) - - xs_modit_lp=modit_xsvector(cnu,indexnu,R,dq,nsigmaD,ngammaL,Sij,nus,ngammaL_grid) - wls_modit = 100000000/nus - - #ref (direct) - d=10000 - ll=mdbCO.nu_lines - xsv_lpf_lp=lpf_xsection(nus,ll,sigmaD,gammaL,Sij,memory_size=30) - - - dif=xs_modit_lp/xsv_lpf_lp-1. - med=np.median(dif) - iju=22940. - ijd=26400. - limu,limd=np.searchsorted(wls_modit[::-1],[iju,ijd]) - std=np.std(dif[::-1][limu:limd]) - - return med,std,R,ijd,iju,wls_modit,xs_modit_lp,xsv_lpf_lp,dif - - -if __name__=="__main__": - import matplotlib - m,std,R,ijd,iju,wls_modit,xs_modit_lp,xsv_lpf_lp,dif=comperr(200000) - m1,std1,R1,ijd1,iju1,wls_modit1,xs_modit_lp1,xsv_lpf_lp1,dif1=comperr(400000) - - print(m,std,R) - print(m1,std1,R1) - - - - - #PLOT - plotfig=True - if plotfig: - matplotlib.rcParams['agg.path.chunksize'] = 100000 - llow=2300.4 - lhigh=2300.7 - tip=2.0 - fig=plt.figure(figsize=(12,3)) - ax=plt.subplot2grid((12, 1), (0, 0),rowspan=8) - plt.plot(wls_modit1,xsv_lpf_lp1,label="Direct",color="C0",alpha=0.3,markersize=3) - plt.plot(wls_modit,xs_modit_lp,color="C1",lw=1,alpha=0.3,label="R="+str(R)) - plt.plot(wls_modit1,xs_modit_lp1,color="C2",lw=1,alpha=0.5,label="R="+str(R1),ls="dashed") - - plt.xlim(ijd,iju) -# plt.xlim(llow*10-tip,lhigh*10+tip) - - plt.ylim(1.1e-35,1.e-17) -# plt.ylim(1.e-27,3.e-17) - plt.yscale("log") - -# plt.xlim(llow*10.0-tip,lhigh*10.0+tip) - plt.legend(loc="upper right") - plt.ylabel(" cross section",fontsize=10) - #plt.text(22986,3.e-21,"$P=10^{-3}$ bar") - plt.xlabel('wavelength [$\AA$]') - - ax=plt.subplot2grid((12, 1), (8, 0),rowspan=4) - plt.plot(wls_modit,np.abs((dif)*100),alpha=0.2,color="C1",label="R="+str(R)) - plt.plot(wls_modit1,np.abs((dif1)*100),alpha=0.5,color="C2",label="R="+str(R1)) - - plt.ylabel("difference (%)",fontsize=10) - plt.xlim(ijd,iju) -# plt.xlim(llow*10-tip,lhigh*10+tip) - - plt.ylim(0.1,10000.0) -# plt.ylim(-10*100*std,10*100*std) - plt.yscale("log") - plt.xlabel('wavelength [$\AA$]') - plt.legend(loc="upper left") - - plt.savefig("fig/comparison_modit.png", bbox_inches="tight", pad_inches=0.0) - plt.savefig("fig/comparison_modit.pdf", bbox_inches="tight", pad_inches=0.0) - plt.show() diff --git a/tests/manual_check/comparison/modit_hitran_CO.py b/tests/manual_check/comparison/modit_hitran_CO.py deleted file mode 100644 index cfd59e07..00000000 --- a/tests/manual_check/comparison/modit_hitran_CO.py +++ /dev/null @@ -1,106 +0,0 @@ -""" MODIT HITRAN CO - -""" -from jax import config -from exojax.spec.hitran import normalized_doppler_sigma -from exojax.spec.set_ditgrid import ditgrid_log_interval -from exojax.spec import api -from exojax.spec.hitran import line_strength, doppler_sigma, gamma_hitran, gamma_natural -from exojax.spec.lpf import auto_xsection as lpf_xsection -from exojax.spec import initspec -from exojax.spec.modit import xsvector as modit_xsvector -from exojax.spec.lpf import xsvector as lpf_xsvector -from exojax.spec import make_numatrix0 -import numpy as np -import tqdm -import jax.numpy as jnp -from jax import vmap -import seaborn as sns -import matplotlib.pyplot as plt -plt.style.use('bmh') - -nus = np.logspace(np.log10(4000), np.log10(4500.0), 3000000, dtype=np.float64) -mdbCO = api.MdbHitran('~/exojax/data/CO/05_hit12.par', nus) - -Mmol = 28.010446441149536 -Tref = 296.0 -Tfix = 1000.0 -Pfix = 1.e-3 - -# USE TIPS partition function -Q296 = np.array([107.25937215917970, 224.38496958496091, 112.61710362499998, - 660.22969049609367, 236.14433662109374, 1382.8672147421873]) -Q1000 = np.array([382.19096582031250, 802.30952197265628, 402.80326733398437, - 2357.1041210937501, 847.84866308593757, 4928.7215078125000]) -qr = Q1000/Q296 - -qt = np.ones_like(mdbCO.isoid, dtype=np.float64) -for idx, iso in enumerate(mdbCO.uniqiso): - mask = mdbCO.isoid == iso - qt[mask] = qr[idx] - -Sij = line_strength(Tfix, mdbCO.logsij0, mdbCO.nu_lines, mdbCO.elower, qt, mdbCO.Tref) -gammaL = gamma_hitran(Pfix, Tfix, Pfix, mdbCO.n_air, - mdbCO.gamma_air, mdbCO.gamma_self) -# + gamma_natural(A) #uncomment if you inclide a natural width -sigmaD = doppler_sigma(mdbCO.nu_lines, Tfix, Mmol) - -cnu, indexnu, R, pmarray = initspec.init_modit(mdbCO.nu_lines, nus) -nsigmaD = normalized_doppler_sigma(Tfix, Mmol, R) -ngammaL = gammaL/(mdbCO.nu_lines/R) -ngammaL_grid = ditgrid_log_interval(ngammaL) - -xs_modit_lp = modit_xsvector( - cnu, indexnu, R, pmarray, nsigmaD, ngammaL, Sij, nus, ngammaL_grid) -wls_modit = 100000000/nus - -#ref (direct) -d = 10000 -ll = mdbCO.nu_lines -xsv_lpf_lp = lpf_xsection(nus, ll, sigmaD, gammaL, Sij, memory_size=30) - -config.update('jax_enable_x64', True) - -xs_modit_lp_f64 = modit_xsvector( - cnu, indexnu, R, pmarray, nsigmaD, ngammaL, Sij, nus, ngammaL_grid) - - -# PLOT -llow = 2300.4 -lhigh = 2300.7 -tip = 20.0 -fig = plt.figure(figsize=(12, 3)) -ax = plt.subplot2grid((12, 1), (0, 0), rowspan=8) -plt.plot(wls_modit, xsv_lpf_lp, label='Direct', - color='C0', markersize=3, alpha=0.3) -plt.plot(wls_modit, xs_modit_lp, lw=1, - color='C1', alpha=1, label='MODIT (F32)') -plt.plot(wls_modit, xs_modit_lp_f64, lw=1, - color='C2', alpha=1, label='MODIT (F64)') - - -plt.ylim(1.1e-28, 1.e-17) -# plt.ylim(1.e-27,3.e-20) -plt.yscale('log') - -plt.xlim(llow*10.0-tip, lhigh*10.0+tip) -plt.legend(loc='upper right') -plt.ylabel(' cross section $(\mathrm{cm}^2)$', fontsize=10) -#plt.text(22986,3.e-21,"$P=10^{-3}$ bar") -plt.xlabel('wavelength [$\AA$]') - -ax = plt.subplot2grid((12, 1), (8, 0), rowspan=4) -plt.plot(wls_modit, np.abs(xs_modit_lp/xsv_lpf_lp-1.)*100, - lw=1, alpha=0.5, color='C1', label='MODIT (F32)') -plt.plot(wls_modit, np.abs(xs_modit_lp_f64/xsv_lpf_lp-1.) * - 100, lw=1, alpha=1, color='C2', label='MODIT (F64)') -plt.yscale('log') -plt.ylabel('difference (%)', fontsize=10) -plt.xlim(llow*10.0-tip, lhigh*10.0+tip) -plt.ylim(0.01, 100.0) -plt.xlabel('wavelength [$\AA$]') -plt.legend(loc='upper left') - -plt.savefig('comparison_modit.png', bbox_inches='tight', pad_inches=0.0) -plt.savefig('comparison_modit.pdf', bbox_inches='tight', pad_inches=0.0) -plt.show()