cloudy: /home66/gary/public_html/cloudy/c08_branch/source/punch

00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
00002  * others.  For conditions of distribution and use see copyright notice in license.txt */
00003 /*punch_opacity punch total opacity in any element, punch opacity command */
00004 #include "cddefines.h"
00005 #include "physconst.h"
00006 #include "iso.h"
00007 #include "rfield.h"
00008 #include "ipoint.h"
00009 #include "continuum.h"
00010 #include "opacity.h"
00011 #include "dense.h"
00012 #include "yield.h"
00013 #include "atmdat.h"
00014 #include "abund.h"
00015 #include "heavy.h"
00016 #include "elementnames.h"
00017 #include "punch.h"
00018 /* print final information about where opacity files are */
00019 STATIC void prtPunOpacSummary(void);
00020 
00021 void punch_opacity(FILE * ioPUN, 
00022   long int ipPun)
00023 {
00024         /* this will be the file name for opacity output */
00025         char chFileName[FILENAME_PATH_LENGTH_2];        
00026 
00027         /* this is its pointer */
00028         FILE *ioFILENAME;
00029 
00030         char chNumbers[31][3] = {"0","1","2","3","4","5","6","7","8","9",
00031                 "10","11","12","13","14","15","16","17","18","19",
00032                 "20","21","22","23","24","25","26","27","28","29","30"};
00033 
00034         long int i, 
00035           ilow, 
00036           ion, 
00037           ipS, 
00038           j, 
00039           nelem;
00040 
00041         double ener, 
00042           ener3;
00043 
00044         DEBUG_ENTRY( "punch_opacity()" );
00045 
00046         /* this flag says to redo the static opacities */
00047         opac.lgRedoStatic = true;
00048 
00049         /* punch total opacity in any element, punch opacity command
00050          * ioPUN is ioPUN unit number, ipPun is pointer within stack of punches */
00051 
00052         if( strcmp(punch.chOpcTyp[ipPun],"TOTL") == 0 )
00053         /* total opacity */
00054         {
00055                 for( j=0; j < rfield.nflux; j++ )
00056                 {
00057                         fprintf( ioPUN, "%12.4e\t%10.2e\t%10.2e\t%10.2e\t%10.2e\t%4.4s\n", 
00058                           AnuUnit(rfield.AnuOrg[j]),
00059                           opac.opacity_abs[j]+opac.OpacStatic[j] + opac.opacity_sct[j], 
00060                           opac.opacity_abs[j]+opac.OpacStatic[j], 
00061                           opac.opacity_sct[j], 
00062                           opac.opacity_sct[j]/MAX2(1e-37,opac.opacity_sct[j]+opac.opacity_abs[j]+opac.OpacStatic[j]), 
00063                           rfield.chContLabel[j] );
00064                 }
00065 
00066                 fprintf( ioPUN, "\n" );
00067         }
00068 
00069         else if( strcmp(punch.chOpcTyp[ipPun],"BREM") == 0 )
00070         /* bremsstrahlung opacity */
00071         {
00072                 for( j=0; j < rfield.nflux; j++ )
00073                 {
00074                         fprintf( ioPUN, "%.5e\t%.3e\n", 
00075                           AnuUnit(rfield.AnuOrg[j]),
00076                           opac.FreeFreeOpacity[j] );
00077                 }
00078 
00079                 fprintf( ioPUN, "\n" );
00080         }
00081 
00082         /* subshell photo cross sections */
00083         else if( strcmp(punch.chOpcTyp[ipPun],"SHEL") == 0 )
00084         {
00085                 nelem = (long)punch.punarg[ipPun][0];
00086                 ion = (long)punch.punarg[ipPun][1];
00087                 ipS = (long)punch.punarg[ipPun][2];
00088                 for( i=opac.ipElement[nelem-1][ion-1][ipS-1][0]-1; 
00089                   i < opac.ipElement[nelem-1][ion-1][ipS-1][1]; i++ )
00090                 {
00091                         fprintf( ioPUN, 
00092                                 "%11.3e %11.3e\n", rfield.anu[i], 
00093                           opac.OpacStack[i-opac.ipElement[nelem-1][ion-1][ipS-1][0]+
00094                           opac.ipElement[nelem-1][ion-1][ipS-1][2]] );
00095                 }
00096         }
00097 
00098         else if( strcmp(punch.chOpcTyp[ipPun],"FINE") == 0 )
00099         {
00100                 /* the fine opacity array */
00101                 realnum sum;
00102                 long int k, nu_hi , nskip;
00103                 if( punch.punarg[ipPun][0] > 0. )
00104                 {
00105                         /* possible lower bounds to energy range - will be zero if not set */
00106                         j = ipFineCont( punch.punarg[ipPun][0] );
00107                 }
00108                 else
00109                 {
00110                         j = 0;
00111                 }
00112 
00113                 /* upper limit */
00114                 if( punch.punarg[ipPun][1]> 0. )
00115                 {
00116                         nu_hi = ipFineCont( punch.punarg[ipPun][1]);
00117                 }
00118                 else
00119                 {
00120                         nu_hi = rfield.nfine;
00121                 }
00122 
00123                 nskip = (long)punch.punarg[ipPun][2];
00124                 ASSERT( nskip > 0 );
00125 
00126                 for( i=j; i<nu_hi; i+=nskip )
00127                 {
00128                         sum = 0;
00129                         for( k=0; k<nskip; ++k )
00130                         {
00131                                 sum += rfield.fine_opac_zone[i+k];
00132                         }
00133                         sum /= nskip;
00134                         if( sum > 0. )
00135                                 fprintf(ioPUN,"%.5e\t%.3e\n",
00136                                         AnuUnit( rfield.fine_anu[i+k/2] ),      sum );
00137                 }
00138         }
00139 
00140         /* figure for hazy */
00141         else if( strcmp(punch.chOpcTyp[ipPun],"FIGU") == 0 )
00142         {
00143                 nelem = 0;
00144                 for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1; i < (iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0] - 1); i++ )
00145                 {
00146                         ener = rfield.anu[i]*0.01356;
00147                         ener3 = 1e24*POW3(ener);
00148                         fprintf( ioPUN, 
00149                                 "%12.4e%12.4e%12.4e%12.4e%12.4e\n", 
00150                           rfield.anu[i], ener, 
00151                           opac.OpacStack[i-iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]+ iso.ipOpac[ipH_LIKE][nelem][ipH1s]]*ener3, 
00152                           0., 
00153                           (opac.opacity_abs[i]+opac.OpacStatic[i])/ dense.gas_phase[ipHYDROGEN]*ener3 );
00154                 }
00155 
00156                 for( i=iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0]-1; i < rfield.nupper; i++ )
00157                 {
00158                         ener = rfield.anu[i]*0.01356;
00159                         ener3 = 1e24*POW3(ener);
00160                         fprintf( ioPUN, 
00161                                 "%12.4e%12.4e%12.4e%12.4e%12.4e\n", 
00162                           rfield.anu[i], 
00163                           ener, 
00164                           opac.OpacStack[i-iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]+ iso.ipOpac[ipH_LIKE][nelem][ipH1s]]*ener3, 
00165                           opac.OpacStack[i-iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0]+ opac.iophe1[0]]*dense.gas_phase[ipHELIUM]/dense.gas_phase[ipHYDROGEN]*ener3, 
00166                           (opac.opacity_abs[i]+opac.OpacStatic[i])/dense.gas_phase[ipHYDROGEN]*ener3 );
00167                 }
00168         }
00169 
00170         /* photoionization data table for AGN */
00171         else if( strcmp(punch.chOpcTyp[ipPun]," AGN") == 0 )
00172         {
00173                 long int 
00174                         ipop,
00175                         nshell,
00176                         nelec;
00177                 char chOutput[100] , chString[100];
00178                 /* now do loop over temp, but add elements */
00179                 for( nelem=ipLITHIUM; nelem<LIMELM; ++nelem )
00180                 {
00181                         /* this list of elements included in the AGN tables is defined in zeroabun.c */
00182                         if( abund.lgAGN[nelem] )
00183                         {
00184                                 for( ion=0; ion<=nelem; ++ion )
00185                                 {
00186                                         /* number of bound electrons */
00187                                         nelec = nelem+1 - ion;
00188 
00189                                         /* print chemical symbol */
00190                                         sprintf(chOutput,"%s", 
00191                                                 elementnames.chElementSym[nelem]);
00192                                         /* some elements have only one letter - this avoids leaving a space */
00193                                         if( chOutput[1]==' ' )
00194                                                 chOutput[1] = chOutput[2];
00195                                         /* now ionization stage */
00196                                         if( ion==0 )
00197                                         {
00198                                                 sprintf(chString,"0 ");
00199                                         }
00200                                         else if( ion==1 )
00201                                         {
00202                                                 sprintf(chString,"+ ");
00203                                         }
00204                                         else
00205                                         {
00206                                                 sprintf(chString,"+%li ",ion);
00207                                         }
00208                                         strcat( chOutput , chString );
00209                                         fprintf(ioPUN,"%s",chOutput );
00210                                         /*fprintf(ioPUN,"\t%.2f\n", Heavy.Valence_IP_Ryd[nelem][ion] );*/
00211 
00212                                         /* now loop over all shells */
00213                                         for( nshell=0; nshell < Heavy.nsShells[nelem][ion]; nshell++ )
00214                                         {
00215                                                 /* shell designation */
00216                                                 fprintf(ioPUN,"\t%s",Heavy.chShell[nshell] );
00217 
00218                                                 /* ionization potential of shell */
00219                                                 fprintf(ioPUN,"\t%.2f" ,
00220                                                         t_ADfA::Inst().ph1(nshell,nelec-1,nelem,0)/EVRYD* 0.9998787);
00221 
00222                                                 /* set lower and upper limits to this range */
00223                                                 ipop = opac.ipElement[nelem][ion][nshell][2];
00224                                                 fprintf(ioPUN,"\t%.2f",opac.OpacStack[ipop-1]/1e-18);
00225                                                 for( i=0; i < t_yield::Inst().nelec_eject(nelem,ion,nshell); ++i )
00226                                                 {
00227                                                         fprintf(ioPUN,"\t%.2f",
00228                                                                 t_yield::Inst().elec_eject_frac(nelem,ion,nshell,i));
00229                                                 }
00230                                                 fprintf(ioPUN,"\n");
00231                                         }
00232 
00233                                 }
00234                         }
00235                 }
00236         }
00237 
00238         /* hydrogen */
00239         else if( strcmp(punch.chOpcTyp[ipPun],"HYDR") == 0 )
00240         {
00241                 nelem = ipHYDROGEN;
00242                 /* zero out the opacity arrays */
00243                 OpacityZero();
00244 
00245                 OpacityAdd1SubshellInduc(iso.ipOpac[ipH_LIKE][nelem][ipH1s],iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s],
00246                   rfield.nupper,1.,0.,'v');
00247                 ilow = Heavy.ipHeavy[nelem][0];
00248 
00249                 /* start filename for output */
00250                 strcpy( chFileName , elementnames.chElementNameShort[0] );
00251 
00252                 /* continue filename with ionization stage */
00253                 strcat( chFileName , chNumbers[1] );
00254 
00255                 /* end it with string ".opc", name will be of form HYDR.opc */
00256                 strcat( chFileName , ".opc" );
00257 
00258                 /* now try to open it for writing */
00259                 ioFILENAME = open_data( chFileName, "w", AS_LOCAL_ONLY );
00260                 for( j=ilow; j <= rfield.nupper; j++ )
00261                 {
00262                         /* photon energy in rydbergs */
00263                         PrintE93( ioFILENAME , rfield.anu[j-1]*EVRYD );
00264                         fprintf( ioFILENAME , "\t");
00265                         /* cross section in megabarns */
00266                         PrintE93( ioFILENAME, (opac.opacity_abs[j-1]+opac.OpacStatic[j-1])/1e-18 );
00267                         fprintf( ioFILENAME , "\n");
00268                 }
00269 
00270                 fclose( ioFILENAME );
00271                 prtPunOpacSummary();
00272                 cdEXIT(EXIT_FAILURE);
00273         }
00274 
00275         /* helium */
00276         else if( strcmp(punch.chOpcTyp[ipPun],"HELI") == 0 )
00277         {
00278                 /* atomic helium first, HELI1.opc */
00279                 nelem = ipHELIUM;
00280                 OpacityZero();
00281                 OpacityAdd1SubshellInduc(opac.iophe1[0],iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0],rfield.nflux,1.,0.,'v');
00282                 ilow = Heavy.ipHeavy[nelem][0];
00283 
00284                 /* start filename for output */
00285                 strcpy( chFileName , elementnames.chElementNameShort[1] );
00286 
00287                 /* continue filename with ionization stage */
00288                 strcat( chFileName , chNumbers[1] );
00289 
00290                 /* end it wil .opc, name will be of form HYDR.opc */
00291                 strcat( chFileName , ".opc" );
00292 
00293                 /* now try to open it for writing */
00294                 ioFILENAME = open_data( chFileName, "w", AS_LOCAL_ONLY );
00295                 for( j=ilow; j <= rfield.nupper; j++ )
00296                 {
00297                         /* photon energy in rydbergs */
00298                         PrintE93( ioFILENAME , rfield.anu[j-1]*EVRYD );
00299                         fprintf( ioFILENAME , "\t");
00300                         /* cross section in megabarns */
00301                         PrintE93( ioFILENAME, (opac.opacity_abs[j-1]+opac.OpacStatic[j-1])/1e-18 );
00302                         fprintf( ioFILENAME , "\n");
00303                 }
00304                 fclose( ioFILENAME );
00305 
00306                 /* now do helium ion, HELI2.opc */
00307                 OpacityZero();
00308                 OpacityAdd1SubshellInduc(iso.ipOpac[ipH_LIKE][1][ipH1s],iso.ipIsoLevNIonCon[ipH_LIKE][1][ipH1s],rfield.nupper,1.,0.,'v');
00309                 ilow = Heavy.ipHeavy[nelem][1];
00310 
00311                 /* start filename for output */
00312                 strcpy( chFileName , elementnames.chElementNameShort[1] );
00313 
00314                 /* continue filename with ionization stage */
00315                 strcat( chFileName , chNumbers[2] );
00316 
00317                 /* end it wil .opc, name will be of form HYDR.opc */
00318                 strcat( chFileName , ".opc" );
00319 
00320                 /* now try to open it for writing */
00321                 ioFILENAME = open_data( chFileName, "w", AS_LOCAL_ONLY );
00322                 for( j=ilow; j <= rfield.nupper; j++ )
00323                 {
00324                         /* photon energy in rydbergs */
00325                         PrintE93( ioFILENAME , rfield.anu[j-1]*EVRYD );
00326                         fprintf( ioFILENAME , "\t");
00327                         /* cross section in megabarns */
00328                         PrintE93( ioFILENAME, (opac.opacity_abs[j-1]+opac.OpacStatic[j-1])/1e-18 );
00329                         fprintf( ioFILENAME , "\n");
00330                 }
00331 
00332                 prtPunOpacSummary();
00333                 fclose( ioFILENAME );
00334                 cdEXIT(EXIT_FAILURE);
00335         }
00336 
00337         else
00338         {
00339                 /* check for hydrogen through zinc, nelem is atomic number on the c scale */
00340                 nelem = -1;
00341                 i = 0;
00342                 while( i < LIMELM )
00343                 {
00344                         if( strcmp(punch.chOpcTyp[ipPun],elementnames.chElementNameShort[i]) == 0 )
00345                         {
00346                                 nelem = i;
00347                                 break;
00348                         }
00349                         ++i;
00350                 }
00351 
00352                 /* nelem is still negative if above loop fell through */
00353                 if( nelem < 0 )
00354                 {
00355                         fprintf( ioQQQ, " Unidentified opacity key=%4.4s\n", 
00356                           punch.chOpcTyp[ipPun] );
00357                         cdEXIT(EXIT_FAILURE);
00358                 }
00359 
00360                 /* pc lint did not pick up the above logice an warned possible negative array index */
00361                 ASSERT( nelem>=0);
00362                 /* generic driving of OpacityAdd1Element */
00363                 StatesElem[ipH_LIKE][nelem][ipH1s].Pop = 1.;
00364                 iso.DepartCoef[ipH_LIKE][nelem][ipH1s] = 0.;
00365                 if( nelem > ipHYDROGEN )
00366                 {
00367                         StatesElem[ipHE_LIKE][nelem][ipH1s].Pop = 1.;
00368                         iso.DepartCoef[ipHE_LIKE][nelem][ipH1s] = 0.;
00369                 }
00370 
00371                 for( ion=0; ion <= nelem; ion++ )
00372                 {
00373                         for( j=0; j < (nelem + 2); j++ )
00374                         {
00375                                 dense.xIonDense[nelem][j] = 0.;
00376                         }
00377 
00378                         dense.xIonDense[nelem][ion] = 1.;
00379 
00380                         OpacityZero();
00381 
00382                         /* generate opacity with standard routine - this is the one
00383                          * called in OpacityAddTotal to make opacities in usual calculations */
00384                         OpacityAdd1Element(nelem);
00385 
00386                         /* start filename for output */
00387                         strcpy( chFileName , elementnames.chElementNameShort[nelem] );
00388 
00389                         /* continue filename with ionization stage */
00390                         strcat( chFileName , chNumbers[ion+1] );
00391 
00392                         /* end it wil .opc, name will be of form HYDR.opc */
00393                         strcat( chFileName , ".opc" );
00394 
00395                         /* now try to open it for writing */
00396                         ioFILENAME = open_data( chFileName, "w", AS_LOCAL_ONLY );
00397 
00398                         ilow = Heavy.ipHeavy[nelem][ion];
00399 
00400                         for( j=ilow-1; j < MIN2(rfield.nflux,continuum.KshellLimit); j++ )
00401                         {
00402                                 /* photon energy in rydbergs */
00403                                 PrintE93( ioFILENAME , rfield.anu[j]*EVRYD );
00404                                 fprintf( ioFILENAME , "\t");
00405 
00406                                 /* cross section in megabarns */
00407                                 PrintE93( ioFILENAME, (opac.opacity_abs[j]+opac.OpacStatic[j])/1e-18 );
00408                                 fprintf( ioFILENAME , "\n");
00409                         }
00410                         /* close this ionization stage */
00411                         fclose( ioFILENAME );
00412                 }
00413 
00414                 prtPunOpacSummary();
00415                 cdEXIT(EXIT_FAILURE);
00416         }
00417 
00418         return;
00419 }
00420 
00421 /* print final information about where opacity files are */
00422 STATIC void prtPunOpacSummary(void)
00423 {
00424         fprintf(ioQQQ,"\n\nThe opacity files have been successfully created.\n");
00425         fprintf(ioQQQ,"The files have names that start with the first 4 characters of the element name.\n");
00426         fprintf(ioQQQ,"There is one file per ion and the number after the element name indicates the ion.\n");
00427         fprintf(ioQQQ,"The energies are in eV and the cross sections in megabarns.\n");
00428         fprintf(ioQQQ,"All end in \".opc\"\n");
00429         fprintf(ioQQQ,"The data only extend to the highest energy in this continuum source.\n");
00430 }
/home66/gary/public_html/cloudy/c08_branch/source/punch_opacity.cpp
