cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
save_opacity.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2022 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 "iso.h"
6 #include "rfield.h"
7 #include "ipoint.h"
8 #include "continuum.h"
9 #include "opacity.h"
10 #include "dense.h"
11 #include "yield.h"
12 #include "atmdat_adfa.h"
13 #include "abund.h"
14 #include "heavy.h"
15 #include "elementnames.h"
16 #include "save.h"
17 #include "freebound.h"
18 #include "species.h"
19 
20 /* print final information about where opacity files are */
21 STATIC void prtPunOpacSummary(void);
22 
23 void save_opacity(FILE * ioPUN,
24  long int ipPun)
25 {
26  /* this will be the file name for opacity output */
27  char chFileName[FILENAME_PATH_LENGTH_2];
28 
29  /* this is its pointer */
30  FILE *ioFILENAME;
31 
32  char chNumbers[31][3] = {"0","1","2","3","4","5","6","7","8","9",
33  "10","11","12","13","14","15","16","17","18","19",
34  "20","21","22","23","24","25","26","27","28","29","30"};
35 
36  long int i,
37  ilow,
38  ion,
39  ipS,
40  j,
41  nelem;
42 
43  double ener,
44  ener3;
45 
46  DEBUG_ENTRY( "save_opacity()" );
47 
48  /* this flag says to redo the static opacities */
49  opac.lgRedoStatic = true;
50 
51  /* save total opacity in any element, save opacity command
52  * ioPUN is ioPUN unit number, ipPun is pointer within stack of punches */
53 
54  if( strcmp(save.chOpcTyp[ipPun],"TOTL") == 0 )
55  /* total opacity */
56  {
57  for( j=0; j < rfield.nflux; j++ )
58  {
59  fprintf( ioPUN, "%12.4e\t%10.2e\t%10.2e\t%10.2e\t%10.2e\t%4.4s\n",
60  AnuUnit(rfield.anu(j)),
63  opac.opacity_sct[j],
65  rfield.chContLabel[j].c_str() );
66  }
67 
68  fprintf( ioPUN, "\n" );
69  }
70 
71  else if( strcmp(save.chOpcTyp[ipPun],"BREM") == 0 )
72  /* bremsstrahlung opacity */
73  {
74  for( j=0; j < rfield.nflux; j++ )
75  {
76  fprintf( ioPUN, "%.5e\t%.3e\t%.3e\n",
77  AnuUnit(rfield.anu(j)),
80  }
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)abs(save.punarg[ipPun][2]);
125  nskip = MAX2( 1, nskip );
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  // ALL key -> report all, even zero
136  if( sum > 0. || save.punarg[ipPun][2]<0)
137  fprintf(ioPUN,"%.5e\t%.3e\n",
138  AnuUnit( rfield.fine_anu[i+k/2] ), sum );
139  }
140  }
141 
142  /* figure for hazy */
143  else if( strcmp(save.chOpcTyp[ipPun],"FIGU") == 0 )
144  {
145  nelem = 0;
146  for( i=iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1; i < (iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon - 1); i++ )
147  {
148  ener = rfield.anu(i)*0.01356;
149  ener3 = 1e24*POW3(ener);
150  fprintf( ioPUN,
151  "%12.4e%12.4e%12.4e%12.4e%12.4e\n",
152  rfield.anu(i), ener,
153  opac.OpacStack[i-iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon+ iso_sp[ipH_LIKE][nelem].fb[ipH1s].ipOpac]*ener3,
154  0.,
156  }
157 
158  for( i=iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon-1; i < rfield.nflux_with_check; i++ )
159  {
160  ener = rfield.anu(i)*0.01356;
161  ener3 = 1e24*POW3(ener);
162  fprintf( ioPUN,
163  "%12.4e%12.4e%12.4e%12.4e%12.4e\n",
164  rfield.anu(i),
165  ener,
166  opac.OpacStack[i-iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon+ iso_sp[ipH_LIKE][nelem].fb[ipH1s].ipOpac]*ener3,
169  }
170  }
171 
172  /* photoionization data table for AGN */
173  else if( strcmp(save.chOpcTyp[ipPun]," AGN") == 0 )
174  {
175  /* now do loop over temp, but add elements */
176  for( nelem=ipLITHIUM; nelem<LIMELM; ++nelem )
177  {
178  /* this list of elements included in the AGN tables is defined in zeroabun.c */
179  if( abund.lgAGN[nelem] )
180  {
181  for( ion=0; ion<=nelem; ++ion )
182  {
183  /* number of bound electrons */
184  long nelec = nelem+1 - ion;
185  fprintf( ioPUN, "%s ", makeChemical( nelem, nelec ).c_str() );
186  /*fprintf(ioPUN,"\t%.2f\n", Heavy.Valence_IP_Ryd[nelem][ion] );*/
187 
188  /* now loop over all shells */
189  for( long nshell=0; nshell < Heavy.nsShells[nelem][ion]; nshell++ )
190  {
191  /* shell designation */
192  fprintf(ioPUN,"\t%s",Heavy.chShell[nshell] );
193 
194  /* ionization potential of shell */
195  fprintf(ioPUN,"\t%.2f" ,
196  t_ADfA::Inst().ph1(nshell,nelec-1,nelem,0)/EVRYD);
197 
198  /* set lower and upper limits to this range */
199  long ipop = opac.ipElement[nelem][ion][nshell][2];
200  fprintf(ioPUN,"\t%.2f",opac.OpacStack[ipop-1]/1e-18);
201  for( i=0; i < t_yield::Inst().nelec_eject(nelem,ion,nshell); ++i )
202  {
203  fprintf(ioPUN,"\t%.2f",
204  t_yield::Inst().elec_eject_frac(nelem,ion,nshell,i));
205  }
206  fprintf(ioPUN,"\n");
207  }
208 
209  }
210  }
211  }
212  }
213 
214  /* hydrogen */
215  else if( strcmp(save.chOpcTyp[ipPun],"HYDR") == 0 )
216  {
217  nelem = ipHYDROGEN;
218  /* zero out the opacity arrays */
219  OpacityZero();
220 
221  OpacityAdd1SubshellInduc(iso_sp[ipH_LIKE][nelem].fb[ipH1s].ipOpac,iso_sp[ipH_LIKE][nelem].fb[ipH1s].ipIsoLevNIonCon,
222  rfield.nflux_with_check,1.,0.,'v');
223  ilow = Heavy.ipHeavy[nelem][0];
224 
225  /* start filename for output */
226  strcpy( chFileName , elementnames.chElementNameShort[0] );
227 
228  /* continue filename with ionization stage */
229  strcat( chFileName , chNumbers[1] );
230 
231  /* end it with string ".opc", name will be of form HYDR.opc */
232  strcat( chFileName , ".opc" );
233 
234  /* now try to open it for writing */
235  ioFILENAME = open_data( chFileName, "w" );
236  for( j=ilow; j <= rfield.nflux_with_check; j++ )
237  {
238  /* photon energy in rydbergs */
239  PrintE93( ioFILENAME , rfield.anu(j-1)*EVRYD );
240  fprintf( ioFILENAME , "\t");
241  /* cross section in megabarns */
242  PrintE93( ioFILENAME, (opac.opacity_abs[j-1]+opac.OpacStatic[j-1])/1e-18 );
243  fprintf( ioFILENAME , "\n");
244  }
245 
246  fclose( ioFILENAME );
249  }
250 
251  /* helium */
252  else if( strcmp(save.chOpcTyp[ipPun],"HELI") == 0 )
253  {
254  /* atomic helium first, HELI1.opc */
255  nelem = ipHELIUM;
256  OpacityZero();
257  OpacityAdd1SubshellInduc(opac.iophe1[0],iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon,rfield.nflux,1.,0.,'v');
258  ilow = Heavy.ipHeavy[nelem][0];
259 
260  /* start filename for output */
261  strcpy( chFileName , elementnames.chElementNameShort[1] );
262 
263  /* continue filename with ionization stage */
264  strcat( chFileName , chNumbers[1] );
265 
266  /* end it wil .opc, name will be of form HYDR.opc */
267  strcat( chFileName , ".opc" );
268 
269  /* now try to open it for writing */
270  ioFILENAME = open_data( chFileName, "w" );
271  for( j=ilow; j <= rfield.nflux_with_check; j++ )
272  {
273  /* photon energy in rydbergs */
274  PrintE93( ioFILENAME , rfield.anu(j-1)*EVRYD );
275  fprintf( ioFILENAME , "\t");
276  /* cross section in megabarns */
277  PrintE93( ioFILENAME, (opac.opacity_abs[j-1]+opac.OpacStatic[j-1])/1e-18 );
278  fprintf( ioFILENAME , "\n");
279  }
280  fclose( ioFILENAME );
281 
282  /* now do helium ion, HELI2.opc */
283  OpacityZero();
284  OpacityAdd1SubshellInduc(iso_sp[ipH_LIKE][1].fb[ipH1s].ipOpac,iso_sp[ipH_LIKE][1].fb[ipH1s].ipIsoLevNIonCon,rfield.nflux_with_check,1.,0.,'v');
285  ilow = Heavy.ipHeavy[nelem][1];
286 
287  /* start filename for output */
288  strcpy( chFileName , elementnames.chElementNameShort[1] );
289 
290  /* continue filename with ionization stage */
291  strcat( chFileName , chNumbers[2] );
292 
293  /* end it wil .opc, name will be of form HYDR.opc */
294  strcat( chFileName , ".opc" );
295 
296  /* now try to open it for writing */
297  ioFILENAME = open_data( chFileName, "w" );
298  for( j=ilow; j <= rfield.nflux_with_check; j++ )
299  {
300  /* photon energy in rydbergs */
301  PrintE93( ioFILENAME , rfield.anu(j-1)*EVRYD );
302  fprintf( ioFILENAME , "\t");
303  /* cross section in megabarns */
304  PrintE93( ioFILENAME, (opac.opacity_abs[j-1]+opac.OpacStatic[j-1])/1e-18 );
305  fprintf( ioFILENAME , "\n");
306  }
307 
309  fclose( ioFILENAME );
311  }
312 
313  else
314  {
315  /* check for hydrogen through zinc, nelem is atomic number on the c scale */
316  nelem = -1;
317  i = 0;
318  while( i < LIMELM )
319  {
320  if( strcmp(save.chOpcTyp[ipPun],elementnames.chElementNameShort[i]) == 0 )
321  {
322  nelem = i;
323  break;
324  }
325  ++i;
326  }
327 
328  /* nelem is still negative if above loop fell through */
329  if( nelem < 0 )
330  {
331  fprintf( ioQQQ, " Unidentified opacity key=%4.4s\n",
332  save.chOpcTyp[ipPun] );
334  }
335 
336  /* pc lint did not pick up the above logice an warned possible negative array index */
337  ASSERT( nelem>=0);
338  /* generic driving of OpacityAdd1Element */
339  iso_sp[ipH_LIKE][nelem].st[ipH1s].Pop() = 1.;
340  iso_sp[ipH_LIKE][nelem].st[ipH1s].DepartCoef() = 0.;
341  if( nelem > ipHYDROGEN )
342  {
343  iso_sp[ipHE_LIKE][nelem].st[ipH1s].Pop() = 1.;
344  iso_sp[ipHE_LIKE][nelem].st[ipH1s].DepartCoef() = 0.;
345  }
346 
347  for( ion=0; ion <= nelem; ion++ )
348  {
349  for( j=0; j < (nelem + 2); j++ )
350  {
351  dense.xIonDense[nelem][j] = 0.;
352  }
353 
354  dense.xIonDense[nelem][ion] = 1.;
355 
356  OpacityZero();
357 
358  /* generate opacity with standard routine - this is the one
359  * called in OpacityAddTotal to make opacities in usual calculations */
360  OpacityAdd1Element(nelem);
361 
362  /* start filename for output */
363  strcpy( chFileName , elementnames.chElementNameShort[nelem] );
364 
365  /* continue filename with ionization stage */
366  strcat( chFileName , chNumbers[ion+1] );
367 
368  /* end it wil .opc, name will be of form HYDR.opc */
369  strcat( chFileName , ".opc" );
370 
371  /* now try to open it for writing */
372  ioFILENAME = open_data( chFileName, "w" );
373 
374  ilow = Heavy.ipHeavy[nelem][ion];
375 
376  for( j=ilow-1; j < MIN2(rfield.nflux,continuum.KshellLimit); j++ )
377  {
378  /* photon energy in rydbergs */
379  PrintE93( ioFILENAME , rfield.anu(j)*EVRYD );
380  fprintf( ioFILENAME , "\t");
381 
382  /* cross section in megabarns */
383  PrintE93( ioFILENAME, (opac.opacity_abs[j]+opac.OpacStatic[j])/1e-18 );
384  fprintf( ioFILENAME , "\n");
385  }
386  /* close this ionization stage */
387  fclose( ioFILENAME );
388  }
389 
392  }
393 
394  return;
395 }
396 
397 /* print final information about where opacity files are */
399 {
400  fprintf(ioQQQ,"\n\nThe opacity files have been successfully created.\n");
401  fprintf(ioQQQ,"The files have names that start with the first 4 characters of the element name.\n");
402  fprintf(ioQQQ,"There is one file per ion and the number after the element name indicates the ion.\n");
403  fprintf(ioQQQ,"The energies are in eV and the cross sections in megabarns.\n");
404  fprintf(ioQQQ,"All end in \".opc\"\n");
405  fprintf(ioQQQ,"The data only extend to the highest energy in this continuum source.\n");
406 }
void OpacityAdd1Element(long int ipZ)
bool lgAGN[LIMELM]
Definition: abund.h:198
realnum * fine_anu
Definition: rfield.h:395
realnum punarg[LIMPUN][3]
Definition: save.h:372
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:765
long int ipElement[LIMELM][LIMELM][7][3]
Definition: opacity.h:223
STATIC void prtPunOpacSummary(void)
void PrintE93(FILE *, double)
Definition: service.cpp:923
double * OpacStack
Definition: opacity.h:164
double * opacity_abs
Definition: opacity.h:104
STATIC long int ipPun
Definition: save_do.cpp:367
qList st
Definition: iso.h:482
const int ipHE_LIKE
Definition: iso.h:65
const int FILENAME_PATH_LENGTH_2
Definition: cddefines.h:296
t_opac opac
Definition: opacity.cpp:5
t_Heavy Heavy
Definition: heavy.cpp:5
double * OpacStatic
Definition: opacity.h:123
vector< string > chContLabel
Definition: rfield.h:213
long int KshellLimit
Definition: continuum.h:98
char chOpcTyp[LIMPUN][5]
Definition: save.h:325
long ipFineCont(double energy_ryd)
FILE * ioQQQ
Definition: cddefines.cpp:7
double * opacity_sct
Definition: opacity.h:107
vector< freeBound > fb
Definition: iso.h:481
#define MIN2(a, b)
Definition: cddefines.h:803
double * eeFreeFreeOpacity
Definition: opacity.h:127
double anu(size_t i) const
Definition: mesh.h:120
t_dense dense
Definition: global.cpp:15
static t_ADfA & Inst()
Definition: cddefines.h:209
t_elementnames elementnames
Definition: elementnames.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
long int nflux_with_check
Definition: rfield.h:49
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
void save_opacity(FILE *io, long int np)
t_abund abund
Definition: abund.cpp:5
long nelec_eject(long n, long i, long ns) const
Definition: yield.h:55
long int nsShells[LIMELM][LIMELM]
Definition: heavy.h:28
const int ipH1s
Definition: iso.h:29
#define STATIC
Definition: cddefines.h:118
t_continuum continuum
Definition: continuum.cpp:6
t_rfield rfield
Definition: rfield.cpp:9
long int iophe1[9]
Definition: opacity.h:223
float realnum
Definition: cddefines.h:124
#define EXIT_FAILURE
Definition: cddefines.h:168
#define cdEXIT(FAIL)
Definition: cddefines.h:482
char chElementNameShort[LIMELM][CHARS_ELEMENT_NAME_SHORT]
Definition: elementnames.h:21
double AnuUnit(realnum energy)
Definition: service.cpp:197
realnum * fine_opac_zone
Definition: rfield.h:391
realnum gas_phase[LIMELM]
Definition: dense.h:76
string makeChemical(long nelem, long ion)
Definition: species.cpp:929
#define ASSERT(exp)
Definition: cddefines.h:613
const int ipH_LIKE
Definition: iso.h:64
const int LIMELM
Definition: cddefines.h:308
T pow2(T a)
Definition: cddefines.h:981
long nfine
Definition: rfield.h:387
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
const int ipHELIUM
Definition: cddefines.h:350
double * FreeFreeOpacity
Definition: opacity.h:126
bool lgRedoStatic
Definition: opacity.h:160
double eden
Definition: dense.h:201
#define MAX2(a, b)
Definition: cddefines.h:824
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
char chShell[7][3]
Definition: heavy.h:31
void OpacityAdd1SubshellInduc(long int ipOpac, long int low, long int ihi, double a, double b, char chStat)
#define POW3
Definition: cddefines.h:986
t_save save
Definition: save.cpp:5
const int ipHYDROGEN
Definition: cddefines.h:349
long int nflux
Definition: rfield.h:46
const int ipLITHIUM
Definition: cddefines.h:351
long int ipHeavy[LIMELM][LIMELM]
Definition: heavy.h:11
void OpacityZero(void)
Definition: opacity_zero.cpp:8
#define EXIT_SUCCESS
Definition: cddefines.h:166