cloudy  trunk
save_opacity.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 /*save_opacity save total opacity in any element, save opacity command */
4 #include "cddefines.h"
5 #include "physconst.h"
6 #include "iso.h"
7 #include "rfield.h"
8 #include "ipoint.h"
9 #include "continuum.h"
10 #include "opacity.h"
11 #include "dense.h"
12 #include "yield.h"
13 #include "atmdat.h"
14 #include "abund.h"
15 #include "heavy.h"
16 #include "elementnames.h"
17 #include "save.h"
18 #include "taulines.h"
19 /* print final information about where opacity files are */
20 STATIC void prtPunOpacSummary(void);
21 
22 void save_opacity(FILE * ioPUN,
23  long int ipPun)
24 {
25  /* this will be the file name for opacity output */
26  char chFileName[FILENAME_PATH_LENGTH_2];
27 
28  /* this is its pointer */
29  FILE *ioFILENAME;
30 
31  char chNumbers[31][3] = {"0","1","2","3","4","5","6","7","8","9",
32  "10","11","12","13","14","15","16","17","18","19",
33  "20","21","22","23","24","25","26","27","28","29","30"};
34 
35  long int i,
36  ilow,
37  ion,
38  ipS,
39  j,
40  nelem;
41 
42  double ener,
43  ener3;
44 
45  DEBUG_ENTRY( "save_opacity()" );
46 
47  /* this flag says to redo the static opacities */
48  opac.lgRedoStatic = true;
49 
50  /* save total opacity in any element, save opacity command
51  * ioPUN is ioPUN unit number, ipPun is pointer within stack of punches */
52 
53  if( strcmp(save.chOpcTyp[ipPun],"TOTL") == 0 )
54  /* total opacity */
55  {
56  for( j=0; j < rfield.nflux; j++ )
57  {
58  fprintf( ioPUN, "%12.4e\t%10.2e\t%10.2e\t%10.2e\t%10.2e\t%4.4s\n",
59  AnuUnit(rfield.AnuOrg[j]),
60  opac.opacity_abs[j]+opac.OpacStatic[j] + opac.opacity_sct[j],
61  opac.opacity_abs[j]+opac.OpacStatic[j],
62  opac.opacity_sct[j],
63  opac.opacity_sct[j]/MAX2(1e-37,opac.opacity_sct[j]+opac.opacity_abs[j]+opac.OpacStatic[j]),
64  rfield.chContLabel[j] );
65  }
66 
67  fprintf( ioPUN, "\n" );
68  }
69 
70  else if( strcmp(save.chOpcTyp[ipPun],"BREM") == 0 )
71  /* bremsstrahlung opacity */
72  {
73  for( j=0; j < rfield.nflux; j++ )
74  {
75  fprintf( ioPUN, "%.5e\t%.3e\n",
76  AnuUnit(rfield.AnuOrg[j]),
77  opac.FreeFreeOpacity[j] );
78  }
79 
80  fprintf( ioPUN, "\n" );
81  }
82 
83  /* subshell photo cross sections */
84  else if( strcmp(save.chOpcTyp[ipPun],"SHEL") == 0 )
85  {
86  nelem = (long)save.punarg[ipPun][0];
87  ion = (long)save.punarg[ipPun][1];
88  ipS = (long)save.punarg[ipPun][2];
89  for( i=opac.ipElement[nelem-1][ion-1][ipS-1][0]-1;
90  i < opac.ipElement[nelem-1][ion-1][ipS-1][1]; i++ )
91  {
92  fprintf( ioPUN,
93  "%11.3e %11.3e\n", rfield.anu[i],
94  opac.OpacStack[i-opac.ipElement[nelem-1][ion-1][ipS-1][0]+
95  opac.ipElement[nelem-1][ion-1][ipS-1][2]] );
96  }
97  }
98 
99  else if( strcmp(save.chOpcTyp[ipPun],"FINE") == 0 )
100  {
101  /* the fine opacity array */
102  realnum sum;
103  long int k, nu_hi , nskip;
104  if( save.punarg[ipPun][0] > 0. )
105  {
106  /* possible lower bounds to energy range - will be zero if not set */
107  j = ipFineCont( save.punarg[ipPun][0] );
108  }
109  else
110  {
111  j = 0;
112  }
113 
114  /* upper limit */
115  if( save.punarg[ipPun][1]> 0. )
116  {
117  nu_hi = ipFineCont( save.punarg[ipPun][1]);
118  }
119  else
120  {
121  nu_hi = rfield.nfine;
122  }
123 
124  nskip = (long)save.punarg[ipPun][2];
125  ASSERT( nskip > 0 );
126 
127  for( i=j; i<nu_hi; i+=nskip )
128  {
129  sum = 0;
130  for( k=0; k<nskip; ++k )
131  {
132  sum += rfield.fine_opac_zone[i+k];
133  }
134  sum /= nskip;
135  if( sum > 0. )
136  fprintf(ioPUN,"%.5e\t%.3e\n",
137  AnuUnit( rfield.fine_anu[i+k/2] ), sum );
138  }
139  }
140 
141  /* figure for hazy */
142  else if( strcmp(save.chOpcTyp[ipPun],"FIGU") == 0 )
143  {
144  nelem = 0;
145  for( i=iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1; i < (iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon - 1); i++ )
146  {
147  ener = rfield.anu[i]*0.01356;
148  ener3 = 1e24*POW3(ener);
149  fprintf( ioPUN,
150  "%12.4e%12.4e%12.4e%12.4e%12.4e\n",
151  rfield.anu[i], ener,
152  opac.OpacStack[i-iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon+ iso_sp[ipH_LIKE][nelem].fb[ipH1s].ipOpac]*ener3,
153  0.,
154  (opac.opacity_abs[i]+opac.OpacStatic[i])/ dense.gas_phase[ipHYDROGEN]*ener3 );
155  }
156 
157  for( i=iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon-1; i < rfield.nupper; i++ )
158  {
159  ener = rfield.anu[i]*0.01356;
160  ener3 = 1e24*POW3(ener);
161  fprintf( ioPUN,
162  "%12.4e%12.4e%12.4e%12.4e%12.4e\n",
163  rfield.anu[i],
164  ener,
165  opac.OpacStack[i-iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon+ iso_sp[ipH_LIKE][nelem].fb[ipH1s].ipOpac]*ener3,
166  opac.OpacStack[i-iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon+ opac.iophe1[0]]*dense.gas_phase[ipHELIUM]/dense.gas_phase[ipHYDROGEN]*ener3,
167  (opac.opacity_abs[i]+opac.OpacStatic[i])/dense.gas_phase[ipHYDROGEN]*ener3 );
168  }
169  }
170 
171  /* photoionization data table for AGN */
172  else if( strcmp(save.chOpcTyp[ipPun]," AGN") == 0 )
173  {
174  long int
175  ipop,
176  nshell,
177  nelec;
178  char chOutput[100] , chString[100];
179  /* now do loop over temp, but add elements */
180  for( nelem=ipLITHIUM; nelem<LIMELM; ++nelem )
181  {
182  /* this list of elements included in the AGN tables is defined in zeroabun.c */
183  if( abund.lgAGN[nelem] )
184  {
185  for( ion=0; ion<=nelem; ++ion )
186  {
187  /* number of bound electrons */
188  nelec = nelem+1 - ion;
189 
190  /* print chemical symbol */
191  sprintf(chOutput,"%s",
192  elementnames.chElementSym[nelem]);
193  /* some elements have only one letter - this avoids leaving a space */
194  if( chOutput[1]==' ' )
195  chOutput[1] = chOutput[2];
196  /* now ionization stage */
197  if( ion==0 )
198  {
199  sprintf(chString,"0 ");
200  }
201  else if( ion==1 )
202  {
203  sprintf(chString,"+ ");
204  }
205  else
206  {
207  sprintf(chString,"+%li ",ion);
208  }
209  strcat( chOutput , chString );
210  fprintf(ioPUN,"%s",chOutput );
211  /*fprintf(ioPUN,"\t%.2f\n", Heavy.Valence_IP_Ryd[nelem][ion] );*/
212 
213  /* now loop over all shells */
214  for( nshell=0; nshell < Heavy.nsShells[nelem][ion]; nshell++ )
215  {
216  /* shell designation */
217  fprintf(ioPUN,"\t%s",Heavy.chShell[nshell] );
218 
219  /* ionization potential of shell */
220  fprintf(ioPUN,"\t%.2f" ,
221  t_ADfA::Inst().ph1(nshell,nelec-1,nelem,0)/EVRYD* 0.9998787);
222 
223  /* set lower and upper limits to this range */
224  ipop = opac.ipElement[nelem][ion][nshell][2];
225  fprintf(ioPUN,"\t%.2f",opac.OpacStack[ipop-1]/1e-18);
226  for( i=0; i < t_yield::Inst().nelec_eject(nelem,ion,nshell); ++i )
227  {
228  fprintf(ioPUN,"\t%.2f",
229  t_yield::Inst().elec_eject_frac(nelem,ion,nshell,i));
230  }
231  fprintf(ioPUN,"\n");
232  }
233 
234  }
235  }
236  }
237  }
238 
239  /* hydrogen */
240  else if( strcmp(save.chOpcTyp[ipPun],"HYDR") == 0 )
241  {
242  nelem = ipHYDROGEN;
243  /* zero out the opacity arrays */
244  OpacityZero();
245 
246  OpacityAdd1SubshellInduc(iso_sp[ipH_LIKE][nelem].fb[ipH1s].ipOpac,iso_sp[ipH_LIKE][nelem].fb[ipH1s].ipIsoLevNIonCon,
247  rfield.nupper,1.,0.,'v');
248  ilow = Heavy.ipHeavy[nelem][0];
249 
250  /* start filename for output */
251  strcpy( chFileName , elementnames.chElementNameShort[0] );
252 
253  /* continue filename with ionization stage */
254  strcat( chFileName , chNumbers[1] );
255 
256  /* end it with string ".opc", name will be of form HYDR.opc */
257  strcat( chFileName , ".opc" );
258 
259  /* now try to open it for writing */
260  ioFILENAME = open_data( chFileName, "w", AS_LOCAL_ONLY );
261  for( j=ilow; j <= rfield.nupper; j++ )
262  {
263  /* photon energy in rydbergs */
264  PrintE93( ioFILENAME , rfield.anu[j-1]*EVRYD );
265  fprintf( ioFILENAME , "\t");
266  /* cross section in megabarns */
267  PrintE93( ioFILENAME, (opac.opacity_abs[j-1]+opac.OpacStatic[j-1])/1e-18 );
268  fprintf( ioFILENAME , "\n");
269  }
270 
271  fclose( ioFILENAME );
272  prtPunOpacSummary();
273  cdEXIT(EXIT_SUCCESS);
274  }
275 
276  /* helium */
277  else if( strcmp(save.chOpcTyp[ipPun],"HELI") == 0 )
278  {
279  /* atomic helium first, HELI1.opc */
280  nelem = ipHELIUM;
281  OpacityZero();
282  OpacityAdd1SubshellInduc(opac.iophe1[0],iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon,rfield.nflux,1.,0.,'v');
283  ilow = Heavy.ipHeavy[nelem][0];
284 
285  /* start filename for output */
286  strcpy( chFileName , elementnames.chElementNameShort[1] );
287 
288  /* continue filename with ionization stage */
289  strcat( chFileName , chNumbers[1] );
290 
291  /* end it wil .opc, name will be of form HYDR.opc */
292  strcat( chFileName , ".opc" );
293 
294  /* now try to open it for writing */
295  ioFILENAME = open_data( chFileName, "w", AS_LOCAL_ONLY );
296  for( j=ilow; j <= rfield.nupper; j++ )
297  {
298  /* photon energy in rydbergs */
299  PrintE93( ioFILENAME , rfield.anu[j-1]*EVRYD );
300  fprintf( ioFILENAME , "\t");
301  /* cross section in megabarns */
302  PrintE93( ioFILENAME, (opac.opacity_abs[j-1]+opac.OpacStatic[j-1])/1e-18 );
303  fprintf( ioFILENAME , "\n");
304  }
305  fclose( ioFILENAME );
306 
307  /* now do helium ion, HELI2.opc */
308  OpacityZero();
309  OpacityAdd1SubshellInduc(iso_sp[ipH_LIKE][1].fb[ipH1s].ipOpac,iso_sp[ipH_LIKE][1].fb[ipH1s].ipIsoLevNIonCon,rfield.nupper,1.,0.,'v');
310  ilow = Heavy.ipHeavy[nelem][1];
311 
312  /* start filename for output */
313  strcpy( chFileName , elementnames.chElementNameShort[1] );
314 
315  /* continue filename with ionization stage */
316  strcat( chFileName , chNumbers[2] );
317 
318  /* end it wil .opc, name will be of form HYDR.opc */
319  strcat( chFileName , ".opc" );
320 
321  /* now try to open it for writing */
322  ioFILENAME = open_data( chFileName, "w", AS_LOCAL_ONLY );
323  for( j=ilow; j <= rfield.nupper; j++ )
324  {
325  /* photon energy in rydbergs */
326  PrintE93( ioFILENAME , rfield.anu[j-1]*EVRYD );
327  fprintf( ioFILENAME , "\t");
328  /* cross section in megabarns */
329  PrintE93( ioFILENAME, (opac.opacity_abs[j-1]+opac.OpacStatic[j-1])/1e-18 );
330  fprintf( ioFILENAME , "\n");
331  }
332 
333  prtPunOpacSummary();
334  fclose( ioFILENAME );
335  cdEXIT(EXIT_SUCCESS);
336  }
337 
338  else
339  {
340  /* check for hydrogen through zinc, nelem is atomic number on the c scale */
341  nelem = -1;
342  i = 0;
343  while( i < LIMELM )
344  {
345  if( strcmp(save.chOpcTyp[ipPun],elementnames.chElementNameShort[i]) == 0 )
346  {
347  nelem = i;
348  break;
349  }
350  ++i;
351  }
352 
353  /* nelem is still negative if above loop fell through */
354  if( nelem < 0 )
355  {
356  fprintf( ioQQQ, " Unidentified opacity key=%4.4s\n",
357  save.chOpcTyp[ipPun] );
358  cdEXIT(EXIT_FAILURE);
359  }
360 
361  /* pc lint did not pick up the above logice an warned possible negative array index */
362  ASSERT( nelem>=0);
363  /* generic driving of OpacityAdd1Element */
364  iso_sp[ipH_LIKE][nelem].st[ipH1s].Pop() = 1.;
365  iso_sp[ipH_LIKE][nelem].fb[ipH1s].DepartCoef = 0.;
366  if( nelem > ipHYDROGEN )
367  {
368  iso_sp[ipHE_LIKE][nelem].st[ipH1s].Pop() = 1.;
369  iso_sp[ipHE_LIKE][nelem].fb[ipH1s].DepartCoef = 0.;
370  }
371 
372  for( ion=0; ion <= nelem; ion++ )
373  {
374  for( j=0; j < (nelem + 2); j++ )
375  {
376  dense.xIonDense[nelem][j] = 0.;
377  }
378 
379  dense.xIonDense[nelem][ion] = 1.;
380 
381  OpacityZero();
382 
383  /* generate opacity with standard routine - this is the one
384  * called in OpacityAddTotal to make opacities in usual calculations */
385  OpacityAdd1Element(nelem);
386 
387  /* start filename for output */
388  strcpy( chFileName , elementnames.chElementNameShort[nelem] );
389 
390  /* continue filename with ionization stage */
391  strcat( chFileName , chNumbers[ion+1] );
392 
393  /* end it wil .opc, name will be of form HYDR.opc */
394  strcat( chFileName , ".opc" );
395 
396  /* now try to open it for writing */
397  ioFILENAME = open_data( chFileName, "w", AS_LOCAL_ONLY );
398 
399  ilow = Heavy.ipHeavy[nelem][ion];
400 
401  for( j=ilow-1; j < MIN2(rfield.nflux,continuum.KshellLimit); j++ )
402  {
403  /* photon energy in rydbergs */
404  PrintE93( ioFILENAME , rfield.anu[j]*EVRYD );
405  fprintf( ioFILENAME , "\t");
406 
407  /* cross section in megabarns */
408  PrintE93( ioFILENAME, (opac.opacity_abs[j]+opac.OpacStatic[j])/1e-18 );
409  fprintf( ioFILENAME , "\n");
410  }
411  /* close this ionization stage */
412  fclose( ioFILENAME );
413  }
414 
415  prtPunOpacSummary();
416  cdEXIT(EXIT_SUCCESS);
417  }
418 
419  return;
420 }
421 
422 /* print final information about where opacity files are */
423 STATIC void prtPunOpacSummary(void)
424 {
425  fprintf(ioQQQ,"\n\nThe opacity files have been successfully created.\n");
426  fprintf(ioQQQ,"The files have names that start with the first 4 characters of the element name.\n");
427  fprintf(ioQQQ,"There is one file per ion and the number after the element name indicates the ion.\n");
428  fprintf(ioQQQ,"The energies are in eV and the cross sections in megabarns.\n");
429  fprintf(ioQQQ,"All end in \".opc\"\n");
430  fprintf(ioQQQ,"The data only extend to the highest energy in this continuum source.\n");
431 }
ipFineCont
long ipFineCont(double energy_ryd)
Definition: cont_ipoint.cpp:226
open_data
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:625
t_continuum::KshellLimit
long int KshellLimit
Definition: continuum.h:122
yield.h
ipHE_LIKE
const int ipHE_LIKE
Definition: iso.h:63
dense
t_dense dense
Definition: dense.cpp:24
elementnames.h
Singleton< t_ADfA >::Inst
static t_ADfA & Inst()
Definition: cddefines.h:175
rfield
t_rfield rfield
Definition: rfield.cpp:8
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
elementnames
t_elementnames elementnames
Definition: elementnames.cpp:5
PrintE93
void PrintE93(FILE *, double)
Definition: service.cpp:838
OpacityZero
void OpacityZero(void)
Definition: opacity_zero.cpp:8
realnum
float realnum
Definition: cddefines.h:103
rfield.h
STATIC
#define STATIC
Definition: cddefines.h:97
ipLITHIUM
const int ipLITHIUM
Definition: cddefines.h:307
abund.h
AS_LOCAL_ONLY
@ AS_LOCAL_ONLY
Definition: cpu.h:208
t_Heavy::nsShells
long int nsShells[LIMELM][LIMELM]
Definition: heavy.h:28
AnuUnit
double AnuUnit(realnum energy)
Definition: service.cpp:173
ipoint.h
t_dense::gas_phase
realnum gas_phase[LIMELM]
Definition: dense.h:71
Heavy
t_Heavy Heavy
Definition: heavy.cpp:5
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
iso.h
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
opac
t_opac opac
Definition: opacity.cpp:5
EXIT_SUCCESS
#define EXIT_SUCCESS
Definition: cddefines.h:138
atmdat.h
MIN2
#define MIN2
Definition: cddefines.h:761
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
t_opac::lgRedoStatic
bool lgRedoStatic
Definition: opacity.h:147
dense.h
prtPunOpacSummary
STATIC void prtPunOpacSummary(void)
Definition: save_opacity.cpp:423
cddefines.h
t_opac::opacity_sct
double * opacity_sct
Definition: opacity.h:98
POW3
#define POW3
Definition: cddefines.h:936
t_elementnames::chElementSym
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
t_opac::OpacStack
double * OpacStack
Definition: opacity.h:151
t_rfield::AnuOrg
double * AnuOrg
Definition: rfield.h:62
t_opac::FreeFreeOpacity
double * FreeFreeOpacity
Definition: opacity.h:117
t_Heavy::chShell
char chShell[7][3]
Definition: heavy.h:31
t_rfield::nflux
long int nflux
Definition: rfield.h:43
heavy.h
t_opac::opacity_abs
double * opacity_abs
Definition: opacity.h:95
t_iso_sp::fb
vector< freeBound > fb
Definition: iso.h:452
MAX2
#define MAX2
Definition: cddefines.h:782
LIMELM
const int LIMELM
Definition: cddefines.h:258
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
t_iso_sp::st
qList st
Definition: iso.h:453
t_yield::nelec_eject
long nelec_eject(long n, long i, long ns) const
Definition: yield.h:55
save.h
t_elementnames::chElementNameShort
char chElementNameShort[LIMELM][CHARS_ELEMENT_NAME_SHORT]
Definition: elementnames.h:21
t_rfield::nupper
long int nupper
Definition: rfield.h:46
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:125
abund
t_abund abund
Definition: abund.cpp:5
t_rfield::fine_opac_zone
realnum * fine_opac_zone
Definition: rfield.h:408
FILENAME_PATH_LENGTH_2
const int FILENAME_PATH_LENGTH_2
Definition: cddefines.h:249
t_rfield::nfine
long nfine
Definition: rfield.h:402
t_opac::OpacStatic
double * OpacStatic
Definition: opacity.h:114
t_rfield::anu
double * anu
Definition: rfield.h:58
t_opac::iophe1
long int iophe1[9]
Definition: opacity.h:233
physconst.h
t_abund::lgAGN
bool lgAGN[LIMELM]
Definition: abund.h:44
OpacityAdd1SubshellInduc
void OpacityAdd1SubshellInduc(long int ipOpac, long int low, long int ihi, double a, double b, char chStat)
Definition: opacity_add1subshell.cpp:65
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:306
taulines.h
t_save::punarg
realnum punarg[LIMPUN][3]
Definition: save.h:254
ipH1s
const int ipH1s
Definition: iso.h:27
continuum
t_continuum continuum
Definition: continuum.cpp:5
opacity.h
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
continuum.h
EVRYD
const UNUSED double EVRYD
Definition: physconst.h:189
t_rfield::fine_anu
realnum * fine_anu
Definition: rfield.h:412
t_save::chOpcTyp
char chOpcTyp[LIMPUN][5]
Definition: save.h:229
t_opac::ipElement
long int ipElement[LIMELM][LIMELM][7][3]
Definition: opacity.h:269
t_Heavy::ipHeavy
long int ipHeavy[LIMELM][LIMELM]
Definition: heavy.h:11
save_opacity
void save_opacity(FILE *ioPUN, long int ipPun)
Definition: save_opacity.cpp:22
OpacityAdd1Element
void OpacityAdd1Element(long int ipZ)
Definition: opacity_add1element.cpp:12
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:62
save
t_save save
Definition: save.cpp:5
t_rfield::chContLabel
char ** chContLabel
Definition: rfield.h:223
Generated by   doxygen 1.8.17