cloudy trunk
Loading...
Searching...
No Matches
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 */
20STATIC void prtPunOpacSummary(void);
21
22void 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 );
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
334 fclose( ioFILENAME );
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] );
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
417 }
418
419 return;
420}
421
422/* print final information about where opacity files are */
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}
t_abund abund
Definition abund.cpp:5
FILE * ioQQQ
Definition cddefines.cpp:7
const int FILENAME_PATH_LENGTH_2
Definition cddefines.h:249
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
#define POW3
Definition cddefines.h:936
#define MIN2
Definition cddefines.h:761
void PrintE93(FILE *, double)
Definition service.cpp:838
const int LIMELM
Definition cddefines.h:258
const int ipLITHIUM
Definition cddefines.h:307
#define EXIT_SUCCESS
Definition cddefines.h:138
#define EXIT_FAILURE
Definition cddefines.h:140
double AnuUnit(realnum energy)
Definition service.cpp:173
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int ipHELIUM
Definition cddefines.h:306
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
static t_ADfA & Inst()
Definition cddefines.h:175
long nelec_eject(long n, long i, long ns) const
Definition yield.h:55
long ipFineCont(double energy_ryd)
t_continuum continuum
Definition continuum.cpp:5
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
@ AS_LOCAL_ONLY
Definition cpu.h:208
t_dense dense
Definition dense.cpp:24
t_elementnames elementnames
t_Heavy Heavy
Definition heavy.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipH1s
Definition iso.h:27
const int ipHE_LIKE
Definition iso.h:63
const int ipH_LIKE
Definition iso.h:62
t_opac opac
Definition opacity.cpp:5
void OpacityZero(void)
void OpacityAdd1SubshellInduc(long int ipOpac, long int low, long int ihi, double a, double b, char chStat)
void OpacityAdd1Element(long int ipZ)
UNUSED const double EVRYD
Definition physconst.h:189
t_rfield rfield
Definition rfield.cpp:8
t_save save
Definition save.cpp:5
STATIC void prtPunOpacSummary(void)
void save_opacity(FILE *ioPUN, long int ipPun)