cloudy
trunk
|
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 }