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

Generated for cloudy by doxygen 1.8.1.1