-
Notifications
You must be signed in to change notification settings - Fork 1
/
model_sed.pro
150 lines (113 loc) · 3.36 KB
/
model_sed.pro
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
PRO MODEL_SED, obj, T, g, metal, ps=ps,eps=eps,over=over,offset=offset,color=color
IF N_params() LT 1 THEN BEGIN
print, 'SYNTAX - MODEL_SED, obj, T, g [metal, nir_scale, /ps, /eps, /over]'
print, 'T AND g can be values or arrays'
RETURN
ENDIF
;calls
;READ_BURROWS, NORM_SPEC, sed_make
; colors_kc
IF N_elements(offset) EQ 0 THEN BEGIN
offset= 0
offset2 = -1
ENDIF ELSE BEGIN
offset2 = 1
ENDELSE
@colors_kc
IF n_elements(color) EQ 0 THEN color=[red,blue,cyan]
obj=strn(obj)
IF keyword_set(ps) OR keyword_set(eps) THEN BEGIN
!p.font=0
set_plot, 'ps'
IF keyword_set(ps) THEN BEGIN
extn='.ps'
device, filename=obj+'_model'+extn, encapsulated=0, /portrait, yoffset=0.5, xoffset=0
ENDIF ELSE BEGIN
extn='.eps'
device, filename=obj+'_model'+extn, encapsulated=1
ENDELSE
device, /helvetica,/color,/inches, xsize=8.0,ysize=10.5
ENDIF
n_g=n_elements(g)
n_T=n_elements(T)
n_m = n_elements(metal)
IF n_g EQ 0 THEN g=[4.5,5,5.5]
IF n_m EQ 0 THEN metal=[0.3,1,3]
n_g=n_elements(g)
n_m = n_elements(metal)
IF ~KEYWORD_SET(over) AND (n_T GE 2 OR n_g GE 2) THEN !p.multi=[0,1,4]
p=0
FOR i= 0, n_T-1 DO BEGIN
n_g=n_elements(g)
n_m = n_elements(metal)
IF T[i] EQ 2100 OR T[i] EQ 2200 THEN BEGIN
g_str='5.0'
n_g=1
metal_s='solar'
n_m=1
ENDIF
FOR j= 0, n_g-1 DO BEGIN
T_str=strn(T[i])
IF T[i] NE 2100 AND T[i] NE 2200 THEN BEGIN
CASE g[j] OF
4.5: g_str='4.5'
5.0: g_str='5.0'
5.5: g_str='5.5'
ENDCASE
ENDIF
n_m = n_elements(metal)
IF T[i] EQ 2000 AND g[j] EQ 5.0 THEN BEGIN
z=where(metal NE 3)
metal2 = metal[z]
n_m = n_elements(metal2)
ENDIF ELSE IF T[i] EQ 2100 OR T[i] EQ 2200 THEN BEGIN
metal2 = 1
n_m=1
ENDIF ELSE BEGIN
metal2=metal
ENDELSE
n_m = n_elements(metal2)
FOR k = 0, n_m-1 DO BEGIN
IF T[i] NE 2100 AND T[i] NE 2200 THEN BEGIN
CASE metal2[k] OF
0.3: metal_s='0.3X'
3: metal_s='3X'
1: metal_s='solar'
ENDCASE
ENDIF
;READ IN MODEL and normalize
model=READ_BURROWS('T'+T_str+'_g'+g_str+'_f100_'+metal_s,/fnu)
model_w=model[0,*]
;Normalize model to K band peak (F_nu)
startw= 2.15 & endw= 2.25
model_f=NORM_SPEC(model_w,model[1,*],startw, endw)
w=where(model_w GT 0.5)
IF KEYWORD_SET(over) THEN BEGIN
;if no offset, draw axes on first plot
IF offset2 EQ -1 AND p EQ 0 THEN make_sed, obj, /fnu,/norm,/noopt
oplot, model_w[w], model_f[w]+offset, color=color[p],thick=1
;XYOUTS, 8, 0.85+offset, 'T= '+T_str+', g= '+ g_str +', metal= '+ metal_s
ENDIF ELSE BEGIN
make_sed, obj, /fnu,/norm,/noopt
oplot, model_w, model_f, color=red,thick=1
XYOUTS, 6, 0.85+offset, 'T= '+T_str+', g= '+ g_str +', metal= '+ metal_s
ENDELSE
p=p+1
IF p MOD 4 EQ 0 AND n_g*n_t*n_m GT 4 AND ~KEYWORD_SET(ps) AND ~KEYWORD_SET(eps) THEN BEGIN
print, ''
Print,'Press any key to continue'
print, ' '
anykey=get_kbrd(1)
ENDIF
ENDFOR ;end loop over metal
ENDFOR ;end loop over g
ENDFOR ;end loop over T
;put data on top of models
IF keyword_set(over) THEN make_sed, obj, /fnu,/norm, offset=offset,/noopt,/grey
IF keyword_set(ps) OR keyword_set(eps) THEN BEGIN
device,/close
set_plot,'x'
PRINT, 'wrote spectrum to file '+ obj+'_model'+extn
ENDIF
!p.multi=0
END