cloudy trunk
Loading...
Searching...
No Matches
cdspec.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/*cdSPEC returns the spectrum needed for Keith Arnaud's XSPEC */
4#include "cddefines.h"
5#include "cddrive.h"
6#include "physconst.h"
7#include "geometry.h"
8#include "radius.h"
9#include "rfield.h"
10#include "opacity.h"
11#include "grid.h"
12
13/*
14 * this routine returns the spectrum needed for Keith Arnaud's XSPEC
15 * X-Ray analysis code. It should be called after cdDrive has successfully
16 * computed a model. the calling routine must ensure that the vectors
17 * have enough space to store the resulting spectrum,
18 * given the bounds and energy resolution
19 */
20
21void cdSPEC(
22 /* option - the type of spectrum to be returned
23 * 1 the incident continuum 4\pi nuJ_nu, , erg cm-2 s-1
24 *
25 * 2 the attenuated incident continuum, same units
26 * 3 the reflected continuum, same units
27 *
28 * 4 diffuse continuous emission outward direction
29 * 5 diffuse continuous emission, reflected
30 *
31 * 6 collisional+recombination lines, outward
32 * 7 collisional+recombination lines, reflected
33 *
34 * 8 pumped lines, outward <= not implemented
35 * 9 pumped lines, reflected <= not implemented
36 *
37 * all lines and continuum emitted by the cloud assume full coverage of
38 * continuum source */
39 int nOption ,
40
41 /* the number of cells + 1*/
42 long int nEnergy ,
43
44 /* the returned spectrum, same size is two energy spectra (see option), returns nEnergy -1 pts */
45 double ReturnedSpectrum[] )
46
47{
48 /* this pointer will bet set to one of the cloudy continuum arrays */
49 realnum *cont ,
50 refac;
51 long int ncell , j;
52
53 /* flag saying whether we will need to free cont at the end */
54 bool lgFREE;
55
56 DEBUG_ENTRY( "cdSPEC()" );
57
58 ASSERT( nEnergy <= rfield.nflux );
59
60 if( nOption == 1 )
61 {
62 /* this is the incident continuum, col 2 of save continuum command */
63 cont = rfield.flux_total_incident[0];
64 lgFREE = false;
65 }
66 else if( nOption == 2 )
67 {
68 /* the attenuated transmitted continuum, no diffuse emission,
69 * col 3 of save continuum command */
70 cont = rfield.flux[0];
71 lgFREE = false;
72 }
73 else if( nOption == 3 )
74 {
75 /* reflected incident continuum, col 6 of save continuum command */
76 lgFREE = false;
77 cont = rfield.ConRefIncid[0];
78 }
79 else if( nOption == 4 )
80 {
81 /* diffuse continuous emission outward direction */
82 cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
83
84 /* need to free the vector once done */
85 lgFREE = true;
86 refac = (realnum)radius.r1r0sq*geometry.covgeo;
87 for( j=0; j<rfield.nflux; ++j)
88 {
89 cont[j] = rfield.ConEmitOut[0][j]*refac;
90 }
91 }
92 else if( nOption == 5 )
93 {
94 /* reflected diffuse continuous emission */
95 cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
96 /* need to free the vector once done */
97 lgFREE = true;
98 refac = (realnum)radius.r1r0sq*geometry.covgeo;
99 for( j=0; j<rfield.nflux; ++j)
100 {
101 cont[j] = rfield.ConEmitReflec[0][j]*refac;
102 }
103 }
104 else if( nOption == 6 )
105 {
106 /* all outward lines, */
107 cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
108 /* need to free the vector once done */
109 lgFREE = true;
110 /* correct back to inner radius */
111 refac = (realnum)radius.r1r0sq*geometry.covgeo;
112 for( j=0; j<rfield.nflux; ++j)
113 {
114 /* units of lines here are to cancel out with tricks applied to continuum cells
115 * when finally extracted below */
116 cont[j] = rfield.outlin[0][j] *rfield.widflx[j]/rfield.anu[j]*refac;
117 }
118 }
119 else if( nOption == 7 )
120 {
121 /* all reflected lines */
122 if( geometry.lgSphere )
123 {
124 refac = 0.;
125 }
126 else
127 {
128 refac = 1.;
129 }
130
131 cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
132 /* need to free the vector once done */
133 lgFREE = true;
134 for( j=0; j<rfield.nflux; ++j)
135 {
136 /* units of lines here are to cancel out with tricks applied to continuum cells
137 * when finally extracted below */
138 cont[j] = rfield.reflin[0][j] *rfield.widflx[j]/rfield.anu[j]*refac;
139 }
140 }
141 else
142 {
143 fprintf(ioQQQ," cdSPEC called with impossible nOption (%i)\n", nOption);
145 }
146
147 /* now generate the continua */
148 for( ncell = 0; ncell < nEnergy-1; ++ncell )
149 {
150 ReturnedSpectrum[ncell] = cont[ncell] * EN1RYD * rfield.anu2[ncell] / rfield.widflx[ncell];
151 }
152
153 /* need to free the vector once done if this flag was set */
154 if( lgFREE )
155 {
156 free(cont);
157 }
158 return;
159}
160
161
162/* returns in units photons/cm^2/s/bin */
164 /* option - the type of spectrum to be returned
165 * 1 the incident continuum 4\pi nuJ_nu, , erg cm-2 s-1
166 *
167 * 2 the attenuated incident continuum, same units
168 * 3 the reflected continuum, same units
169 *
170 * 4 diffuse emission, lines + continuum, outward
171 * 5 diffuse emission, lines + continuum, reflected
172 *
173 * 6 diffuse continuous emission outward direction
174 * 7 diffuse continuous emission, reflected
175 *
176 * 8 total transmitted, lines and continuum
177 * 9 total reflected, lines and continuum
178 *
179 *10 exp(-tau) to the illuminated face
180 *
181 * all lines and continuum emitted by the cloud assume full coverage of
182 * continuum source */
183 int nOption ,
184
185 /* the number of cells */
186 long int nEnergy,
187
188 /* the index of the lowest and highest energy bounds to use. */
189 long ipLoEnergy,
190 long ipHiEnergy,
191
192 /* the returned spectrum, same size is two energy spectra (see option), returns nEnergy -1 pts */
193 realnum ReturnedSpectrum[] )
194
195{
196 realnum refac;
197
198 DEBUG_ENTRY( "cdSPEC2()" );
199
200 ASSERT( ipLoEnergy >= 0 );
201 ASSERT( ipLoEnergy < ipHiEnergy );
202 ASSERT( ipHiEnergy < rfield.nupper );
203 ASSERT( nEnergy == (ipHiEnergy-ipLoEnergy+1) );
204 ASSERT( nEnergy >= 2 );
205
206 ASSERT( nOption <= NUM_OUTPUT_TYPES );
207
208 const realnum *trans_coef_total = rfield.getCoarseTransCoef();
209
210 for( long i = 0; i < nEnergy; i++ )
211 {
212 long j = ipLoEnergy + i;
213
214 if( j >= rfield.nflux )
215 {
216 ReturnedSpectrum[i] = SMALLFLOAT;
217 continue;
218 }
219
220 if( nOption == 0 )
221 {
222 /* the attenuated incident continuum */
223 realnum flxatt = rfield.flux[0][j]*
224 (realnum)radius.r1r0sq * trans_coef_total[j];
225
226 /* the outward emitted continuum */
227 realnum conem = (rfield.ConEmitOut[0][j] + rfield.outlin[0][j])*
228 (realnum)radius.r1r0sq*geometry.covgeo;
229
230 /* the reflected continuum */
231 realnum flxref = rfield.ConRefIncid[0][j] + rfield.ConEmitReflec[0][j] +
232 rfield.reflin[0][j];
233
234 ReturnedSpectrum[i] = flxatt + conem + flxref;
235 }
236 else if( nOption == 1 )
237 {
238 /* this is the incident continuum, col 2 of save continuum command */
239 ReturnedSpectrum[i] = rfield.flux_total_incident[0][j];
240 }
241 else if( nOption == 2 )
242 {
243 /* the attenuated transmitted continuum, no diffuse emission,
244 * col 3 of save continuum command */
245 ReturnedSpectrum[i] = rfield.flux[0][j]*
246 (realnum)radius.r1r0sq * trans_coef_total[j];
247 }
248 else if( nOption == 3 )
249 {
250 /* reflected incident continuum, col 6 of save continuum command */
251 ReturnedSpectrum[i] = rfield.ConRefIncid[0][j];
252 }
253 else if( nOption == 4 )
254 {
255 /* all outward diffuse emission */
256 /* correct back to inner radius */
257 refac = (realnum)radius.r1r0sq*geometry.covgeo;
258 ReturnedSpectrum[i] = (rfield.outlin[0][j]+rfield.ConEmitOut[0][j])*refac;
259 }
260 else if( nOption == 5 )
261 {
262 /* all reflected diffuse emission */
263 if( geometry.lgSphere )
264 {
265 refac = 0.;
266 }
267 else
268 {
269 refac = 1.;
270 }
271
272 ReturnedSpectrum[i] = (rfield.reflin[0][j]+rfield.ConEmitReflec[0][j])*refac;
273 }
274 else if( nOption == 6 )
275 {
276 /* all outward line emission */
277 /* correct back to inner radius */
278 refac = (realnum)radius.r1r0sq*geometry.covgeo;
279 ReturnedSpectrum[i] = rfield.outlin[0][j]*refac;
280 }
281 else if( nOption == 7 )
282 {
283 /* all reflected line emission */
284 if( geometry.lgSphere )
285 {
286 refac = 0.;
287 }
288 else
289 {
290 refac = 1.;
291 }
292
293 ReturnedSpectrum[i] = rfield.reflin[0][j]*refac;
294 }
295 else if( nOption == 8 )
296 {
297 /* total transmitted continuum */
298 /* correct back to inner radius */
299 refac = (realnum)radius.r1r0sq*geometry.covgeo;
300 ReturnedSpectrum[i] = (rfield.ConEmitOut[0][j]+ rfield.outlin[0][j])*refac
301 + rfield.flux[0][j]*(realnum)radius.r1r0sq*trans_coef_total[j];
302 }
303 else if( nOption == 9 )
304 {
305 /* total reflected continuum */
306 ReturnedSpectrum[i] = rfield.ConRefIncid[0][j] + rfield.ConEmitReflec[0][j] +
307 rfield.reflin[0][j];
308 }
309 else if( nOption == 10 )
310 {
311 /* this is exp(-tau) */
312 /* This is the TOTAL attenuation in both the continuum and lines.
313 * Jon Miller discovered that the line attenuation was missing in 07.02 */
314 ReturnedSpectrum[i] = opac.ExpmTau[j]*trans_coef_total[j];
315 }
316 else
317 {
318 fprintf(ioQQQ," cdSPEC called with impossible nOption (%i)\n", nOption);
320 }
321
322 ASSERT( ReturnedSpectrum[i] >=0.f );
323 }
324
325 return;
326}
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
#define MALLOC(exp)
Definition cddefines.h:501
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void cdSPEC2(int nOption, long int nEnergy, long ipLoEnergy, long ipHiEnergy, realnum ReturnedSpectrum[])
Definition cdspec.cpp:163
void cdSPEC(int nOption, long int nEnergy, double ReturnedSpectrum[])
Definition cdspec.cpp:21
const realnum SMALLFLOAT
Definition cpu.h:191
t_geometry geometry
Definition geometry.cpp:5
const int NUM_OUTPUT_TYPES
Definition grid.h:21
t_opac opac
Definition opacity.cpp:5
UNUSED const double EN1RYD
Definition physconst.h:179
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8