Science/cloudy-agn
1
0
Fork 0
mirror of https://asciireactor.com/otho/cloudy-agn.git synced 2025-10-31 21:28:06 +00:00
Code Issues Packages Projects Releases Wiki Activity

working on magdziarz sed program

Browse Source
This commit is contained in:
caes 2016-06-01 18:51:48 -04:00
parent a6bea59974
commit 25d52ba30b
2 changed files with 27 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
scripts/xmgrace/SED_EFe.par → scripts/xmgrace/SED_rel.par
View File

@ -84,7 +84,7 @@ timestamp color 1
timestamp rot 0 timestamp rot 0
timestamp font 0 timestamp font 0
timestamp char size 1.000000 timestamp char size 1.000000
timestamp def "Wed Jun 1 02:40:17 2016" timestamp def "Wed Jun 1 12:59:43 2016"
r0 off r0 off
link r0 to g0 link r0 to g0
r0 type above r0 type above
@ -131,7 +131,7 @@ g0 fixedpoint xy 0.000000, 0.000000
g0 fixedpoint format general general g0 fixedpoint format general general
g0 fixedpoint prec 6, 6 g0 fixedpoint prec 6, 6
with g0 with g0
world 0.001, 1e-40, 1000, 10000 world 0.001, 0.01, 100, 4
stack world 0, 0, 0, 0 stack world 0, 0, 0, 0
znorm 1 znorm 1
view 0.150000, 0.150000, 1.150000, 0.850000 view 0.150000, 0.150000, 1.150000, 0.850000
@ -205,7 +205,7 @@ with g0
yaxis bar color 1 yaxis bar color 1
yaxis bar linestyle 1 yaxis bar linestyle 1
yaxis bar linewidth 1.0 yaxis bar linewidth 1.0
yaxis label "nuFnu (erg / cm^2 / s )" yaxis label "Relative Value"
yaxis label layout para yaxis label layout para
yaxis label place auto yaxis label place auto
yaxis label char size 1.000000 yaxis label char size 1.000000

40
src/sed/ngc5548_magdziarz.py
View File

@ -1,10 +1,7 @@
""" A program to recreate the average spectrum fit for ngc 5548 in Magdziarz et al 1998 """
#parents, babies = (1, 1)
#while babies < 100:
# #print 'This generation has {0} babies'.format(babies)
# print('This generation has',babies,'babies.')
# parents, babies = (babies, parents + babies)
import math import math
import scipy
import numpy
PI=3.14159265358979323846; PI=3.14159265358979323846;
PLANCK_CONST=4.135668e-15; # in eV * s PLANCK_CONST=4.135668e-15; # in eV * s
@ -13,8 +10,8 @@ RYDBERG_CONST=1.0973731568539e7; # in 1 / m
RYDBERG_UNIT_EV=13.60569252; # in eV RYDBERG_UNIT_EV=13.60569252; # in eV
RYDBERG_UNIT_ANGSTROM=1e10/RYDBERG_CONST; # in A RYDBERG_UNIT_ANGSTROM=1e10/RYDBERG_CONST; # in A
CONT_MIN_ENERGY_keV = 1e-5; # eV CONT_MIN_ENERGY_keV = 1e-3;
CONT_MAX_ENERGY_keV = 1e2; # eV CONT_MAX_ENERGY_keV = 1e2;
CONT_MIN_X = math.log10(CONT_MIN_ENERGY_keV); CONT_MIN_X = math.log10(CONT_MIN_ENERGY_keV);
CONT_MAX_X = math.log10(CONT_MAX_ENERGY_keV); CONT_MAX_X = math.log10(CONT_MAX_ENERGY_keV);
CONT_WIDTH_X = CONT_MAX_X - CONT_MIN_X; CONT_WIDTH_X = CONT_MAX_X - CONT_MIN_X;
@ -41,8 +38,8 @@ IN_EV_2500A = 12398.41929/2500;
# Soft Excess # Soft Excess
α_SE = 1.1 # Quoted consistent with Korista 1995 and Marshall 1997 α_SE_sec2_1 = 1.1 # Quoted consistent with Korista 1995 and Marshall 1997
kT_SE = .56 kT_SE_sec2_1 = .56
# Comtonization fitted to ROSAT data # Comtonization fitted to ROSAT data
@ -51,12 +48,17 @@ kT_SE = .56
# OSSE data fit # OSSE data fit
α_HC = 0.86 α_HC = 0.86
R = 0.96 R = 0.96
E_cutoff_HC = 120 # keV, phase 1 E_cutoff_HC = .120 # keV, phase 1
F_HC = .38 # keV cm⁻² s⁻¹ F_HC = .38 # keV cm⁻² s⁻¹
# E_cutoff_HC = 118 # keV, phase 3 # E_cutoff_HC = 118 # keV, phase 3
# F_HC = .61 # keV cm⁻² s⁻¹ # F_HC = .61 # keV cm⁻² s⁻¹
# Section 3.3 values
kT_SE_sec3_2 = .270 # keV
α_SE_sec3_2 = 1.13
kT_HC_sec3_2 = 55 # keV
α_HC_sec3_2 = .76
@ -78,7 +80,7 @@ def histogram_table(n):
output.append(value) output.append(value)
# Add a final point at 100 KeV # Add a final point at 100 KeV
hν = 1e5; hν = 1e2;
value = sed(hν); value = sed(hν);
output.append((hν,value)) output.append((hν,value))
return output; return output;
@ -86,6 +88,9 @@ def histogram_table(n):
def sed(hν): def sed(hν):
magnitude=0.0; magnitude=0.0;
magnitude += powlaw_cutoff(hν,α_HC,E_cutoff_HC,1) # OSSE data fit magnitude += powlaw_cutoff(hν,α_HC,E_cutoff_HC,1) # OSSE data fit
# magnitude += powlaw_cutoff(hν,α_SE_sec2_1,kT_SE_sec2_1,1)
# magnitude += powlaw_cutoff(hν,α_SE_sec3_2,kT_SE_sec3_2,1)
#magnitude += compt_approx(hν,-1.3,.345,.0034,1)
if magnitude < CONT_MIN_VAL: return CONT_MIN_VAL if magnitude < CONT_MIN_VAL: return CONT_MIN_VAL
# magnitude = CONT_MIN_VAL; # magnitude = CONT_MIN_VAL;
return magnitude; return magnitude;
@ -94,13 +99,16 @@ def powlaw_cutoff(hν,α,E_cutoff,norm):
low_cutoff = .1 low_cutoff = .1
resultant = norm resultant = norm
resultant *= math.exp(-hν/E_cutoff) resultant *= math.exp(-hν/E_cutoff)
resultant *= math.exp(-low_cutoff/hν) #resultant *= math.exp(-low_cutoff/hν)
resultant *= math.pow(hν,1+α) resultant *= math.pow(hν,1+α)
return resultant return resultant
# for t in (22.6, 25.8, 27.3, 29.8):+P B def compt_approx(hν,α,kT_keV,cutoff_keV,norm):
# print(t, ": ", fahrenheit(t)) magnitude = math.pow(hν,(1+α))
magnitude *= math.exp(-(hν/kT_keV))
magnitude *= math.exp(-(cutoff_keV/hν))
magnitude *= norm
return magnitude
test_table = histogram_table(500) test_table = histogram_table(500)