cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cdspec.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 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 /*Spec_cont called by cdSPEC to generate actual spectrum */
5 #include "cddefines.h"
6 #include "cddrive.h"
7 #include "physconst.h"
8 #include "geometry.h"
9 #include "radius.h"
10 #include "rfield.h"
11 #include "opacity.h"
12 #include "grid.h"
13 
14 /*
15  * this routine returns the spectrum needed for Keith Arnaud's XSPEC
16  * X-Ray analysis code. It should be called after cdDrive has successfully
17  * computed a model. the calling routine must ensure that the vectors
18  * have enough space to store the resulting spectrum,
19  * given the bounds and energy resolution
20  */
21 
22 /* array index within energy array in Cloudy grids - this has file scope */
23 static long int iplo , iphi;
24 
25 /* energies of lower and upper bound of XSPEC cell we are considering */
26 static double Elo , Ehi;
27 
28 /* interpolate on cloudy's predicted spectrum to get results on XSPEC grid */
29 STATIC double Spec_cont( realnum cont[] );
30 
31 void cdSPEC(
32  /* option - the type of spectrum to be returned
33  * 1 the incident continuum 4\pi nuJ_nu, , erg cm-2 s-1
34  *
35  * 2 the attenuated incident continuum, same units
36  * 3 the reflected continuum, same units
37  *
38  * 4 diffuse continuous emission outward direction
39  * 5 diffuse continuous emission, reflected
40  *
41  * 6 collisional+recombination lines, outward
42  * 7 collisional+recombination lines, reflected
43  *
44  * 8 pumped lines, outward <= not implemented
45  * 9 pumped lines, reflected <= not implemented
46  *
47  * all lines and continuum emitted by the cloud assume full coverage of
48  * continuum source */
49  int nOption ,
50 
51  /* the energy of the lower edge of each cell
52  * (in Ryd to be consistent with all the rest of Cloudy) */
53  double EnergyLow[] ,
54 
55  /* the number of cells + 1*/
56  long int nEnergy ,
57 
58  /* the returned spectrum, same size is two energy spectra (see option), returns nEnergy -1 pts */
59  double ReturnedSpectrum[] )
60 
61 {
62  /* this pointer will bet set to one of the cloudy continuum arrays */
63  realnum *cont ,
64  refac;
65  long int ncell , j;
66 
67  /* flag saying whether we will need to free cont at the end */
68  bool lgFREE;
69 
70  DEBUG_ENTRY( "cdSPEC()" );
71 
72  if( nOption == 1 )
73  {
74  /* this is the incident continuum, col 2 of punch continuum command */
76  lgFREE = false;
77  }
78  else if( nOption == 2 )
79  {
80  /* the attenuated transmitted continuum, no diffuse emission,
81  * col 3 of punch continuum command */
82  cont = rfield.flux;
83  lgFREE = false;
84  }
85  else if( nOption == 3 )
86  {
87  /* reflected incident continuum, col 6 of punch continuum command */
88  lgFREE = false;
89  cont = rfield.ConRefIncid;
90  }
91  else if( nOption == 4 )
92  {
93  /* diffuse continuous emission outward direction */
94  cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
95 
96  /* need to free the vector once done */
97  lgFREE = true;
99  for( j=0; j<rfield.nflux; ++j)
100  {
101  cont[j] = rfield.ConEmitOut[j]*refac;
102  }
103  }
104  else if( nOption == 5 )
105  {
106  /* reflected diffuse continuous emission */
107  cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
108  /* need to free the vector once done */
109  lgFREE = true;
111  for( j=0; j<rfield.nflux; ++j)
112  {
113  cont[j] = rfield.ConEmitReflec[j]*refac;
114  }
115  }
116  else if( nOption == 6 )
117  {
118  /* all outward lines, */
119  cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
120  /* need to free the vector once done */
121  lgFREE = true;
122  /* correct back to inner radius */
124  for( j=0; j<rfield.nflux; ++j)
125  {
126  /* units of lines here are to cancel out with tricks applied to continuum cells
127  * when finally extracted below */
128  cont[j] = rfield.outlin[j] *rfield.widflx[j]/rfield.anu[j]*refac;
129  }
130  }
131  else if( nOption == 7 )
132  {
133  /* all reflected lines */
134  if( geometry.lgSphere )
135  {
136  refac = 0.;
137  }
138  else
139  {
140  refac = 1.;
141  }
142 
143  cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
144  /* need to free the vector once done */
145  lgFREE = true;
146  for( j=0; j<rfield.nflux; ++j)
147  {
148  /* units of lines here are to cancel out with tricks applied to continuum cells
149  * when finally extracted below */
150  cont[j] = rfield.reflin[j] *rfield.widflx[j]/rfield.anu[j]*refac;
151  }
152  }
153  else
154  {
155  fprintf(ioQQQ," cdSPEC called with impossible nOption (%i)\n", nOption);
156  cdEXIT(EXIT_FAILURE);
157  }
158 
159  /* this is array index used in Spec_cont */
160  iplo = 0;
161  iphi = 0;
162  /* now generate the continua */
163  for( ncell = 0; ncell < nEnergy-1; ++ncell )
164  {
165  /* Spec_cont will find the spectrum between these two energies */
166  Elo = EnergyLow[ncell];
167  Ehi = EnergyLow[ncell+1];
168  ReturnedSpectrum[ncell] = Spec_cont( cont );
169  }
170 
171  /* need to free the vector once done if this flag was set */
172  if( lgFREE )
173  {
174  free(cont);
175  }
176  return;
177 }
178 
179 /* interpolate on cloudy's predicted spectrum to get results on XSPEC grid */
180 STATIC double Spec_cont( realnum cont[] )
181 {
182  double sum;
183  long int i;
184 
185  DEBUG_ENTRY( "Spec_cont()" );
186 
187  /* iplo and iphi are set to zero before this is called to generate a full spectrum,
188  * since energies are in increasing order, after the first call it should only
189  * be necessary to increase ip above where it now is */
190  /* return zero if continuum does not go this high */
191  if( iplo >= rfield.nflux )
192  return 0.;
193 
194  /* now skip out to where continuum starts, often we will already be there
195  * if this is a successive call */
196  /* iplo should be index where Elo is greater than anu-widflx/1. */
197  while( (iplo < rfield.nflux) &&
198  !(Elo <= (rfield.AnuOrg[iplo+1]-rfield.widflx[iplo+1]/2.) &&
199  (Elo >= (rfield.AnuOrg[iplo]-rfield.widflx[iplo]/2.) )) )
200  ++iplo;
201 
202  iphi = iplo;
203  /* iphi should be index where Ehi is greater than anu-widflx/1. */
204  while( (iphi < rfield.nflux) &&
205  !(Ehi <= (rfield.AnuOrg[iphi+1]-rfield.widflx[iphi+1]/2.) &&
206  (Ehi >= (rfield.AnuOrg[iphi]-rfield.widflx[iphi]/2.) )) )
207  ++iphi;
208 
209  sum = 0.;
210  if( iphi-iplo>1 )
211  {
212  /* this branch when more than two cells are involved -
213  * begin by summing intermediate cells */
214  for( i=iplo+1; i<iphi; ++i )
215  {
216  sum += cont[i] * rfield.anu2[i];
217  }
218  }
219  ASSERT( sum>=0.);
220 
221  /* at this point sum will often be zero, and we need to add on flux from
222  * the iplo and iphi cells - there are two possibilities - iplo and iphi can
223  * be the same cell, or different cells */
224  if( iplo == iphi )
225  {
226  /* we want only a piece of the ip cell */
227  sum = cont[iplo]/rfield.widflx[iplo]*(Ehi-Elo) * rfield.anu2[iplo];
228  ASSERT( sum>=0.);
229  }
230  else
231  {
232  /* we want to add on a fraction of the low and high energy cells */
233  double piece;
234 
235  /* the low energy cell */
236  piece = cont[iplo]/rfield.widflx[iplo]*( (rfield.AnuOrg[iplo+1]-rfield.widflx[iplo+1]/2.) - Elo ) *
237  rfield.anu2[iplo];
238  ASSERT( piece >= 0.);
239  sum += piece;
240 
241  /* now the high energy cell */
242  piece = cont[iphi]/rfield.widflx[iphi]*(Ehi - (rfield.AnuOrg[iphi]-rfield.widflx[iphi]/2.) ) *
243  rfield.anu2[iphi];
244  ASSERT( piece >= 0.);
245  sum += piece;
246  ASSERT( sum>=0.);
247  }
248 
249  /* now divide by total energy width */
250  sum *= EN1RYD / (Ehi - Elo );
251 
252  ASSERT( sum >=0. );
253  return sum;
254 
255 # if 0
256  /* at this point Elo is less than the upper edge energy of the ip cell */
257  /* >>chng 01 jul 23, change logic, use this branch if desired width smaller than
258  * intrinsic cloudy width */
259  /*if( rfield.anu[ip]+rfield.widflx[ip]/2. > Ehi )*/
260  if( (Ehi-Elo) < rfield.widflx[ip] )
261  {
262  /* this is where we want the flux evaluated */
263  double Emiddle = (Elo + Ehi)/2.;
264  double f1 , f2 ,finterp;
265  ASSERT( (Elo >= (rfield.anu[ip]+rfield.widflx[ip]/2.
266  /* this is case where requested cell is within a single Cloudy cell -
267  * do linear interpolation */
268  /* >>chng 01 jul 23, change to linear interpolation */
269  /*return( cont[ip] /rfield.widflx[ip] *rfield.anu2[ip]*EN1RYD );*/
270  f1 = cont[ip-1] /rfield.widflx[ip-1] *rfield.anu2[ip-1]*EN1RYD;
271  f2 = cont[ip] /rfield.widflx[ip] *rfield.anu2[ip]*EN1RYD;
272 
273  finterp = f1 + (f2-f1) /rfield.widflx[ip] *
274  fabs(Emiddle-(rfield.anu[ip]-rfield.widflx[ip]/2.));
275  return( finterp );
276  }
277  else
278  {
279  /* this is case where we will add on more than one Cloudy cell */
280  sum = cont[ip] * rfield.anu2[ip]* EN1RYD;
281  /* energy of upper edge of Cloudy cell */
282  EcHi = ( rfield.anu[ip] + rfield.widflx[ip]/2. + rfield.anu[ip+1] - rfield.widflx[ip+1]/2.)/2.;
283  /* mutiply this by energy width of this first cell that counts to final integral */
284  sum *= (EcHi - Elo )/rfield.widflx[ip];
285  }
286 
287  /* increment index since we will now work on next cell */
288  ++ip;
289 
290  /* now add up total cells, this loop for cloudy cells totally within desired bins,
291  * we will divide by energy width later, so factor of widflx is missing here */
292  while( Ehi > rfield.anu[ip]+rfield.widflx[ip]/2. )
293  {
294  sum += cont[ip] * rfield.anu2[ip]*EN1RYD;
295  ++ip;
296  }
297 
298  /* energy of upper edge of Cloudy cell */
299  EcLo = ( rfield.anu[ip] - rfield.widflx[ip]/2. + rfield.anu[ip-1] + rfield.widflx[ip-1]/2.)/2.;
300 
301  /* now finish up with last cell, whose upper bound is higher than desired bin */
302  sum += cont[ip] * rfield.anu2[ip] * EN1RYD *
303  (Ehi - EcLo )/rfield.widflx[ip];
304 
305  /* now divide by total energy width */
306  sum /= (Ehi - Elo );
307 
308  return(sum);
309 # endif
310 }
311 
312 /* returns in units photons/cm^2/s/bin */
313 void cdSPEC2(
314  /* option - the type of spectrum to be returned
315  * 1 the incident continuum 4\pi nuJ_nu, , erg cm-2 s-1
316  *
317  * 2 the attenuated incident continuum, same units
318  * 3 the reflected continuum, same units
319  *
320  * 4 diffuse continuous emission outward direction
321  * 5 diffuse continuous emission, reflected
322  *
323  * 6 collisional+recombination lines, outward
324  * 7 collisional+recombination lines, reflected
325  *
326  * 8 total transmitted, lines and continuum
327  * 9 total reflected, lines and continuum
328  *
329  *10 exp(-tau) to the illuminated face
330  *
331  * all lines and continuum emitted by the cloud assume full coverage of
332  * continuum source */
333  int nOption ,
334 
335  /* the number of cells + 1*/
336  long int nEnergy ,
337 
338  /* the returned spectrum, same size is two energy spectra (see option), returns nEnergy -1 pts */
339  realnum ReturnedSpectrum[] )
340 
341 {
342  /* this pointer will bet set to one of the cloudy continuum arrays */
343  realnum *cont ,
344  refac;
345  long int ncell , j;
346 
347  /* flag saying whether we will need to free cont at the end */
348  bool lgFREE;
349 
350  DEBUG_ENTRY( "cdSPEC2()" );
351 
352  ASSERT( nOption <= NUM_OUTPUT_TYPES );
353 
354  if( nOption == 1 )
355  {
356  /* this is the incident continuum, col 2 of punch continuum command */
358  lgFREE = false;
359  }
360  else if( nOption == 2 )
361  {
362  /* the attenuated transmitted continuum, no diffuse emission,
363  * col 3 of punch continuum command */
364  cont = rfield.flux;
365  lgFREE = false;
366  }
367  else if( nOption == 3 )
368  {
369  /* reflected incident continuum, col 6 of punch continuum command */
370  lgFREE = false;
371  cont = rfield.ConRefIncid;
372  }
373  else if( nOption == 4 )
374  {
375  /* diffuse continuous emission outward direction */
376  cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
377  /* need to free the vector once done */
378  lgFREE = true;
380  for( j=0; j<rfield.nflux; ++j)
381  {
382  cont[j] = rfield.ConEmitOut[j]*refac;
383  }
384  }
385  else if( nOption == 5 )
386  {
387  /* reflected diffuse continuous emission */
388  cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
389  /* need to free the vector once done */
390  lgFREE = true;
391  for( j=0; j<rfield.nflux; ++j)
392  {
393  cont[j] = rfield.ConEmitReflec[j];
394  }
395  }
396  else if( nOption == 6 )
397  {
398  /* all outward lines, */
399  cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
400  /* need to free the vector once done */
401  lgFREE = true;
402  /* correct back to inner radius */
404  for( j=0; j<rfield.nflux; ++j)
405  {
406  cont[j] = rfield.outlin[j]*refac;
407  }
408  }
409  else if( nOption == 7 )
410  {
411  /* all reflected lines */
412  if( geometry.lgSphere )
413  {
414  refac = 0.;
415  }
416  else
417  {
418  refac = 1.;
419  }
420 
421  cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
422  /* need to free the vector once done */
423  lgFREE = true;
424  for( j=0; j<rfield.nflux; ++j)
425  {
426  /* units of lines here are to cancel out with tricks applied to continuum cells
427  * when finally extracted below */
428  cont[j] = rfield.reflin[j]*refac;
429  }
430  }
431  else if( nOption == 8 )
432  {
433  /* total transmitted continuum */
434  cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
435  /* need to free the vector once done */
436  lgFREE = true;
437  /* correct back to inner radius */
439  for( j=0; j<rfield.nflux; ++j)
440  {
441  /* units of lines here are to cancel out with tricks applied to continuum cells
442  * when finally extracted below */
443  cont[j] = (rfield.ConEmitOut[j]+ rfield.outlin[j])*refac
445  }
446  }
447  else if( nOption == 9 )
448  {
449  /* total reflected continuum */
450  cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
451  /* need to free the vector once done */
452  lgFREE = true;
453  for( j=0; j<rfield.nflux; ++j)
454  {
455  /* units of lines here are to cancel out with tricks applied to continuum cells
456  * when finally extracted below */
457  cont[j] = ((rfield.ConRefIncid[j]+rfield.ConEmitReflec[j]) +
458  rfield.reflin[j]);
459  }
460  }
461  else if( nOption == 10 )
462  {
463  /* this is exp(-tau) */
464  cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
465  /* need to free the vector once done */
466  lgFREE = true;
467  for( j=0; j<nEnergy; ++j)
468  {
469  /* This is the TOTAL attenuation in both the continuum and lines.
470  * Jon Miller discovered that the line attenuation was missing in 07.02 */
471  cont[j] = opac.ExpmTau[j]*rfield.trans_coef_total[j];
472  }
473  }
474  else
475  {
476  fprintf(ioQQQ," cdSPEC called with impossible nOption (%i)\n", nOption);
477  cdEXIT(EXIT_FAILURE);
478  }
479 
480  /* now generate the continua */
481  for( ncell = 0; ncell < nEnergy; ++ncell )
482  {
483  if( ncell >= rfield.nflux )
484  {
485  ReturnedSpectrum[ncell] = SMALLFLOAT;
486  }
487  else
488  {
489  ReturnedSpectrum[ncell] = cont[ncell];
490  }
491  ASSERT( ReturnedSpectrum[ncell] >=0.f );
492  }
493 
494  /* need to free the vector once done if this flag was set */
495  if( lgFREE )
496  {
497  free(cont);
498  }
499  return;
500 }

Generated for cloudy by doxygen 1.8.1.1