00001
00002
00003
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
00019 static void prtPunOpacSummary(void);
00020
00021 void punch_opacity(FILE * ioPUN,
00022 long int ipPun)
00023 {
00024
00025 char chFileName[FILENAME_PATH_LENGTH_2];
00026
00027
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
00047 opac.lgRedoStatic = true;
00048
00049
00050
00051
00052 if( strcmp(punch.chOpcTyp[ipPun],"TOTL") == 0 )
00053
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
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
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
00101 float sum;
00102 long int k, nu_hi , nskip;
00103 if( punch.punarg[ipPun][0] > 0. )
00104 {
00105
00106 j = ipFineCont( punch.punarg[ipPun][0] );
00107 }
00108 else
00109 {
00110 j = 0;
00111 }
00112
00113
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
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
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
00179 for( nelem=ipLITHIUM; nelem<LIMELM; ++nelem )
00180 {
00181
00182 if( abund.lgAGN[nelem] )
00183 {
00184 for( ion=0; ion<=nelem; ++ion )
00185 {
00186
00187 nelec = nelem+1 - ion;
00188
00189
00190 sprintf(chOutput,"%s",
00191 elementnames.chElementSym[nelem]);
00192
00193 if( chOutput[1]==' ' )
00194 chOutput[1] = chOutput[2];
00195
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
00211
00212
00213 for( nshell=0; nshell < Heavy.nsShells[nelem][ion]; nshell++ )
00214 {
00215
00216 fprintf(ioPUN,"\t%s",Heavy.chShell[nshell] );
00217
00218
00219 fprintf(ioPUN,"\t%.2f" ,
00220 t_ADfA::Inst().ph1(nshell,nelec-1,nelem,0)/EVRYD* 0.9998787);
00221
00222
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
00239 else if( strcmp(punch.chOpcTyp[ipPun],"HYDR") == 0 )
00240 {
00241 nelem = ipHYDROGEN;
00242
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
00250 strcpy( chFileName , elementnames.chElementNameShort[0] );
00251
00252
00253 strcat( chFileName , chNumbers[1] );
00254
00255
00256 strcat( chFileName , ".opc" );
00257
00258
00259 if( (ioFILENAME = fopen( chFileName , "w" )) ==NULL)
00260 {
00261 fprintf( ioQQQ," could not open opacity file %s for writing.\n",chFileName);
00262 puts( "[Stop in punopac]" );
00263 cdEXIT(EXIT_FAILURE);
00264 }
00265 ASSERT( ioFILENAME != NULL );
00266 for( j=ilow; j <= rfield.nupper; j++ )
00267 {
00268
00269 PrintE93( ioFILENAME , rfield.anu[j-1]*EVRYD );
00270 fprintf( ioFILENAME , "\t");
00271
00272 PrintE93( ioFILENAME, (opac.opacity_abs[j-1]+opac.OpacStatic[j-1])/1e-18 );
00273 fprintf( ioFILENAME , "\n");
00274 }
00275
00276 fclose( ioFILENAME );
00277 prtPunOpacSummary();
00278 puts( "[Stop in punopac]" );
00279 cdEXIT(EXIT_FAILURE);
00280 }
00281
00282
00283 else if( strcmp(punch.chOpcTyp[ipPun],"HELI") == 0 )
00284 {
00285
00286 nelem = ipHELIUM;
00287 OpacityZero();
00288 OpacityAdd1SubshellInduc(opac.iophe1[0],iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0],rfield.nflux,1.,0.,'v');
00289 ilow = Heavy.ipHeavy[nelem][0];
00290
00291
00292 strcpy( chFileName , elementnames.chElementNameShort[1] );
00293
00294
00295 strcat( chFileName , chNumbers[1] );
00296
00297
00298 strcat( chFileName , ".opc" );
00299
00300
00301 if( NULL==(ioFILENAME = fopen( chFileName , "w" )) )
00302 {
00303 fprintf( ioQQQ," could not open opacity file %s for writing.\n",chFileName);
00304 puts( "[Stop in punopac]" );
00305 cdEXIT(EXIT_FAILURE);
00306 }
00307
00308 ASSERT( ioFILENAME != NULL );
00309 for( j=ilow; j <= rfield.nupper; j++ )
00310 {
00311
00312 PrintE93( ioFILENAME , rfield.anu[j-1]*EVRYD );
00313 fprintf( ioFILENAME , "\t");
00314
00315 PrintE93( ioFILENAME, (opac.opacity_abs[j-1]+opac.OpacStatic[j-1])/1e-18 );
00316 fprintf( ioFILENAME , "\n");
00317 }
00318 fclose( ioFILENAME );
00319
00320
00321 OpacityZero();
00322 OpacityAdd1SubshellInduc(iso.ipOpac[ipH_LIKE][1][ipH1s],iso.ipIsoLevNIonCon[ipH_LIKE][1][ipH1s],rfield.nupper,1.,0.,'v');
00323 ilow = Heavy.ipHeavy[nelem][1];
00324
00325
00326 strcpy( chFileName , elementnames.chElementNameShort[1] );
00327
00328
00329 strcat( chFileName , chNumbers[2] );
00330
00331
00332 strcat( chFileName , ".opc" );
00333
00334
00335 if( NULL==(ioFILENAME = fopen( chFileName , "w" )) )
00336 {
00337 fprintf( ioQQQ," could not open opacity file %s for writing.\n",chFileName);
00338 puts( "[Stop in punopac]" );
00339 cdEXIT(EXIT_FAILURE);
00340 }
00341
00342 ASSERT( ioFILENAME != NULL );
00343 for( j=ilow; j <= rfield.nupper; j++ )
00344 {
00345
00346 PrintE93( ioFILENAME , rfield.anu[j-1]*EVRYD );
00347 fprintf( ioFILENAME , "\t");
00348
00349 PrintE93( ioFILENAME, (opac.opacity_abs[j-1]+opac.OpacStatic[j-1])/1e-18 );
00350 fprintf( ioFILENAME , "\n");
00351 }
00352
00353 prtPunOpacSummary();
00354 fclose( ioFILENAME );
00355 puts( "[Stop in punopac]" );
00356 cdEXIT(EXIT_FAILURE);
00357 }
00358
00359 else
00360 {
00361
00362 nelem = -1;
00363 i = 0;
00364 while( i < LIMELM )
00365 {
00366 if( strcmp(punch.chOpcTyp[ipPun],elementnames.chElementNameShort[i]) == 0 )
00367 {
00368 nelem = i;
00369 break;
00370 }
00371 ++i;
00372 }
00373
00374
00375 if( nelem < 0 )
00376 {
00377 fprintf( ioQQQ, " Unidentified opacity key=%4.4s\n",
00378 punch.chOpcTyp[ipPun] );
00379 puts( "[Stop in punopac]" );
00380 cdEXIT(EXIT_FAILURE);
00381 }
00382
00383
00384 ASSERT( nelem>=0);
00385
00386 iso.Pop2Ion[ipH_LIKE][nelem][ipH1s] = 1.;
00387 iso.DepartCoef[ipH_LIKE][nelem][ipH1s] = 0.;
00388 if( nelem > ipHYDROGEN )
00389 {
00390 iso.Pop2Ion[ipHE_LIKE][nelem][ipH1s] = 1.;
00391 iso.DepartCoef[ipHE_LIKE][nelem][ipH1s] = 0.;
00392 }
00393
00394 for( ion=0; ion <= nelem; ion++ )
00395 {
00396 for( j=0; j < (nelem + 2); j++ )
00397 {
00398 dense.xIonDense[nelem][j] = 0.;
00399 }
00400
00401 dense.xIonDense[nelem][ion] = 1.;
00402
00403 OpacityZero();
00404
00405
00406
00407 OpacityAdd1Element(nelem);
00408
00409
00410 strcpy( chFileName , elementnames.chElementNameShort[nelem] );
00411
00412
00413 strcat( chFileName , chNumbers[ion+1] );
00414
00415
00416 strcat( chFileName , ".opc" );
00417
00418
00419 if( NULL==(ioFILENAME = fopen( chFileName , "w" )) )
00420 {
00421 fprintf( ioQQQ," could not open opacity file %s for writing.\n",chFileName);
00422 puts( "[Stop in punopac]" );
00423 cdEXIT(EXIT_FAILURE);
00424 }
00425
00426 ilow = Heavy.ipHeavy[nelem][ion];
00427
00428 ASSERT( ioFILENAME != NULL );
00429 for( j=ilow-1; j < MIN2(rfield.nflux,continuum.KshellLimit); j++ )
00430 {
00431
00432 PrintE93( ioFILENAME , rfield.anu[j]*EVRYD );
00433 fprintf( ioFILENAME , "\t");
00434
00435
00436 PrintE93( ioFILENAME, (opac.opacity_abs[j]+opac.OpacStatic[j])/1e-18 );
00437 fprintf( ioFILENAME , "\n");
00438 }
00439
00440 fclose( ioFILENAME );
00441 }
00442
00443 prtPunOpacSummary();
00444 puts( "[Stop in punopac]" );
00445 cdEXIT(EXIT_FAILURE);
00446 }
00447
00448
00449 DEBUG_EXIT( "punch_opacity()" );
00450 return;
00451 }
00452
00453
00454 static void prtPunOpacSummary(void)
00455 {
00456 fprintf(ioQQQ,"\n\nThe opacity files have been successfully created.\n");
00457 fprintf(ioQQQ,"The files have names that start with the first 4 characters of the element name.\n");
00458 fprintf(ioQQQ,"There is one file per ion and the number after the element name indicates the ion.\n");
00459 fprintf(ioQQQ,"The energies are in eV and the cross sections in megabarns.\n");
00460 fprintf(ioQQQ,"All end in \".opc\"\n");
00461 fprintf(ioQQQ,"The data only extend to the highest energy in this continuum source.\n");
00462 }