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 "save.h"
00018 #include "taulines.h"
00019
00020 STATIC void prtPunOpacSummary(void);
00021
00022 void save_opacity(FILE * ioPUN,
00023 long int ipPun)
00024 {
00025
00026 char chFileName[FILENAME_PATH_LENGTH_2];
00027
00028
00029 FILE *ioFILENAME;
00030
00031 char chNumbers[31][3] = {"0","1","2","3","4","5","6","7","8","9",
00032 "10","11","12","13","14","15","16","17","18","19",
00033 "20","21","22","23","24","25","26","27","28","29","30"};
00034
00035 long int i,
00036 ilow,
00037 ion,
00038 ipS,
00039 j,
00040 nelem;
00041
00042 double ener,
00043 ener3;
00044
00045 DEBUG_ENTRY( "save_opacity()" );
00046
00047
00048 opac.lgRedoStatic = true;
00049
00050
00051
00052
00053 if( strcmp(save.chOpcTyp[ipPun],"TOTL") == 0 )
00054
00055 {
00056 for( j=0; j < rfield.nflux; j++ )
00057 {
00058 fprintf( ioPUN, "%12.4e\t%10.2e\t%10.2e\t%10.2e\t%10.2e\t%4.4s\n",
00059 AnuUnit(rfield.AnuOrg[j]),
00060 opac.opacity_abs[j]+opac.OpacStatic[j] + opac.opacity_sct[j],
00061 opac.opacity_abs[j]+opac.OpacStatic[j],
00062 opac.opacity_sct[j],
00063 opac.opacity_sct[j]/MAX2(1e-37,opac.opacity_sct[j]+opac.opacity_abs[j]+opac.OpacStatic[j]),
00064 rfield.chContLabel[j] );
00065 }
00066
00067 fprintf( ioPUN, "\n" );
00068 }
00069
00070 else if( strcmp(save.chOpcTyp[ipPun],"BREM") == 0 )
00071
00072 {
00073 for( j=0; j < rfield.nflux; j++ )
00074 {
00075 fprintf( ioPUN, "%.5e\t%.3e\n",
00076 AnuUnit(rfield.AnuOrg[j]),
00077 opac.FreeFreeOpacity[j] );
00078 }
00079
00080 fprintf( ioPUN, "\n" );
00081 }
00082
00083
00084 else if( strcmp(save.chOpcTyp[ipPun],"SHEL") == 0 )
00085 {
00086 nelem = (long)save.punarg[ipPun][0];
00087 ion = (long)save.punarg[ipPun][1];
00088 ipS = (long)save.punarg[ipPun][2];
00089 for( i=opac.ipElement[nelem-1][ion-1][ipS-1][0]-1;
00090 i < opac.ipElement[nelem-1][ion-1][ipS-1][1]; i++ )
00091 {
00092 fprintf( ioPUN,
00093 "%11.3e %11.3e\n", rfield.anu[i],
00094 opac.OpacStack[i-opac.ipElement[nelem-1][ion-1][ipS-1][0]+
00095 opac.ipElement[nelem-1][ion-1][ipS-1][2]] );
00096 }
00097 }
00098
00099 else if( strcmp(save.chOpcTyp[ipPun],"FINE") == 0 )
00100 {
00101
00102 realnum sum;
00103 long int k, nu_hi , nskip;
00104 if( save.punarg[ipPun][0] > 0. )
00105 {
00106
00107 j = ipFineCont( save.punarg[ipPun][0] );
00108 }
00109 else
00110 {
00111 j = 0;
00112 }
00113
00114
00115 if( save.punarg[ipPun][1]> 0. )
00116 {
00117 nu_hi = ipFineCont( save.punarg[ipPun][1]);
00118 }
00119 else
00120 {
00121 nu_hi = rfield.nfine;
00122 }
00123
00124 nskip = (long)save.punarg[ipPun][2];
00125 ASSERT( nskip > 0 );
00126
00127 for( i=j; i<nu_hi; i+=nskip )
00128 {
00129 sum = 0;
00130 for( k=0; k<nskip; ++k )
00131 {
00132 sum += rfield.fine_opac_zone[i+k];
00133 }
00134 sum /= nskip;
00135 if( sum > 0. )
00136 fprintf(ioPUN,"%.5e\t%.3e\n",
00137 AnuUnit( rfield.fine_anu[i+k/2] ), sum );
00138 }
00139 }
00140
00141
00142 else if( strcmp(save.chOpcTyp[ipPun],"FIGU") == 0 )
00143 {
00144 nelem = 0;
00145 for( i=iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1; i < (iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon - 1); i++ )
00146 {
00147 ener = rfield.anu[i]*0.01356;
00148 ener3 = 1e24*POW3(ener);
00149 fprintf( ioPUN,
00150 "%12.4e%12.4e%12.4e%12.4e%12.4e\n",
00151 rfield.anu[i], ener,
00152 opac.OpacStack[i-iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon+ iso_sp[ipH_LIKE][nelem].fb[ipH1s].ipOpac]*ener3,
00153 0.,
00154 (opac.opacity_abs[i]+opac.OpacStatic[i])/ dense.gas_phase[ipHYDROGEN]*ener3 );
00155 }
00156
00157 for( i=iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon-1; i < rfield.nupper; i++ )
00158 {
00159 ener = rfield.anu[i]*0.01356;
00160 ener3 = 1e24*POW3(ener);
00161 fprintf( ioPUN,
00162 "%12.4e%12.4e%12.4e%12.4e%12.4e\n",
00163 rfield.anu[i],
00164 ener,
00165 opac.OpacStack[i-iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon+ iso_sp[ipH_LIKE][nelem].fb[ipH1s].ipOpac]*ener3,
00166 opac.OpacStack[i-iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon+ opac.iophe1[0]]*dense.gas_phase[ipHELIUM]/dense.gas_phase[ipHYDROGEN]*ener3,
00167 (opac.opacity_abs[i]+opac.OpacStatic[i])/dense.gas_phase[ipHYDROGEN]*ener3 );
00168 }
00169 }
00170
00171
00172 else if( strcmp(save.chOpcTyp[ipPun]," AGN") == 0 )
00173 {
00174 long int
00175 ipop,
00176 nshell,
00177 nelec;
00178 char chOutput[100] , chString[100];
00179
00180 for( nelem=ipLITHIUM; nelem<LIMELM; ++nelem )
00181 {
00182
00183 if( abund.lgAGN[nelem] )
00184 {
00185 for( ion=0; ion<=nelem; ++ion )
00186 {
00187
00188 nelec = nelem+1 - ion;
00189
00190
00191 sprintf(chOutput,"%s",
00192 elementnames.chElementSym[nelem]);
00193
00194 if( chOutput[1]==' ' )
00195 chOutput[1] = chOutput[2];
00196
00197 if( ion==0 )
00198 {
00199 sprintf(chString,"0 ");
00200 }
00201 else if( ion==1 )
00202 {
00203 sprintf(chString,"+ ");
00204 }
00205 else
00206 {
00207 sprintf(chString,"+%li ",ion);
00208 }
00209 strcat( chOutput , chString );
00210 fprintf(ioPUN,"%s",chOutput );
00211
00212
00213
00214 for( nshell=0; nshell < Heavy.nsShells[nelem][ion]; nshell++ )
00215 {
00216
00217 fprintf(ioPUN,"\t%s",Heavy.chShell[nshell] );
00218
00219
00220 fprintf(ioPUN,"\t%.2f" ,
00221 t_ADfA::Inst().ph1(nshell,nelec-1,nelem,0)/EVRYD* 0.9998787);
00222
00223
00224 ipop = opac.ipElement[nelem][ion][nshell][2];
00225 fprintf(ioPUN,"\t%.2f",opac.OpacStack[ipop-1]/1e-18);
00226 for( i=0; i < t_yield::Inst().nelec_eject(nelem,ion,nshell); ++i )
00227 {
00228 fprintf(ioPUN,"\t%.2f",
00229 t_yield::Inst().elec_eject_frac(nelem,ion,nshell,i));
00230 }
00231 fprintf(ioPUN,"\n");
00232 }
00233
00234 }
00235 }
00236 }
00237 }
00238
00239
00240 else if( strcmp(save.chOpcTyp[ipPun],"HYDR") == 0 )
00241 {
00242 nelem = ipHYDROGEN;
00243
00244 OpacityZero();
00245
00246 OpacityAdd1SubshellInduc(iso_sp[ipH_LIKE][nelem].fb[ipH1s].ipOpac,iso_sp[ipH_LIKE][nelem].fb[ipH1s].ipIsoLevNIonCon,
00247 rfield.nupper,1.,0.,'v');
00248 ilow = Heavy.ipHeavy[nelem][0];
00249
00250
00251 strcpy( chFileName , elementnames.chElementNameShort[0] );
00252
00253
00254 strcat( chFileName , chNumbers[1] );
00255
00256
00257 strcat( chFileName , ".opc" );
00258
00259
00260 ioFILENAME = open_data( chFileName, "w", AS_LOCAL_ONLY );
00261 for( j=ilow; j <= rfield.nupper; j++ )
00262 {
00263
00264 PrintE93( ioFILENAME , rfield.anu[j-1]*EVRYD );
00265 fprintf( ioFILENAME , "\t");
00266
00267 PrintE93( ioFILENAME, (opac.opacity_abs[j-1]+opac.OpacStatic[j-1])/1e-18 );
00268 fprintf( ioFILENAME , "\n");
00269 }
00270
00271 fclose( ioFILENAME );
00272 prtPunOpacSummary();
00273 cdEXIT(EXIT_SUCCESS);
00274 }
00275
00276
00277 else if( strcmp(save.chOpcTyp[ipPun],"HELI") == 0 )
00278 {
00279
00280 nelem = ipHELIUM;
00281 OpacityZero();
00282 OpacityAdd1SubshellInduc(opac.iophe1[0],iso_sp[ipHE_LIKE][ipHELIUM].fb[0].ipIsoLevNIonCon,rfield.nflux,1.,0.,'v');
00283 ilow = Heavy.ipHeavy[nelem][0];
00284
00285
00286 strcpy( chFileName , elementnames.chElementNameShort[1] );
00287
00288
00289 strcat( chFileName , chNumbers[1] );
00290
00291
00292 strcat( chFileName , ".opc" );
00293
00294
00295 ioFILENAME = open_data( chFileName, "w", AS_LOCAL_ONLY );
00296 for( j=ilow; j <= rfield.nupper; j++ )
00297 {
00298
00299 PrintE93( ioFILENAME , rfield.anu[j-1]*EVRYD );
00300 fprintf( ioFILENAME , "\t");
00301
00302 PrintE93( ioFILENAME, (opac.opacity_abs[j-1]+opac.OpacStatic[j-1])/1e-18 );
00303 fprintf( ioFILENAME , "\n");
00304 }
00305 fclose( ioFILENAME );
00306
00307
00308 OpacityZero();
00309 OpacityAdd1SubshellInduc(iso_sp[ipH_LIKE][1].fb[ipH1s].ipOpac,iso_sp[ipH_LIKE][1].fb[ipH1s].ipIsoLevNIonCon,rfield.nupper,1.,0.,'v');
00310 ilow = Heavy.ipHeavy[nelem][1];
00311
00312
00313 strcpy( chFileName , elementnames.chElementNameShort[1] );
00314
00315
00316 strcat( chFileName , chNumbers[2] );
00317
00318
00319 strcat( chFileName , ".opc" );
00320
00321
00322 ioFILENAME = open_data( chFileName, "w", AS_LOCAL_ONLY );
00323 for( j=ilow; j <= rfield.nupper; j++ )
00324 {
00325
00326 PrintE93( ioFILENAME , rfield.anu[j-1]*EVRYD );
00327 fprintf( ioFILENAME , "\t");
00328
00329 PrintE93( ioFILENAME, (opac.opacity_abs[j-1]+opac.OpacStatic[j-1])/1e-18 );
00330 fprintf( ioFILENAME , "\n");
00331 }
00332
00333 prtPunOpacSummary();
00334 fclose( ioFILENAME );
00335 cdEXIT(EXIT_SUCCESS);
00336 }
00337
00338 else
00339 {
00340
00341 nelem = -1;
00342 i = 0;
00343 while( i < LIMELM )
00344 {
00345 if( strcmp(save.chOpcTyp[ipPun],elementnames.chElementNameShort[i]) == 0 )
00346 {
00347 nelem = i;
00348 break;
00349 }
00350 ++i;
00351 }
00352
00353
00354 if( nelem < 0 )
00355 {
00356 fprintf( ioQQQ, " Unidentified opacity key=%4.4s\n",
00357 save.chOpcTyp[ipPun] );
00358 cdEXIT(EXIT_FAILURE);
00359 }
00360
00361
00362 ASSERT( nelem>=0);
00363
00364 iso_sp[ipH_LIKE][nelem].st[ipH1s].Pop() = 1.;
00365 iso_sp[ipH_LIKE][nelem].fb[ipH1s].DepartCoef = 0.;
00366 if( nelem > ipHYDROGEN )
00367 {
00368 iso_sp[ipHE_LIKE][nelem].st[ipH1s].Pop() = 1.;
00369 iso_sp[ipHE_LIKE][nelem].fb[ipH1s].DepartCoef = 0.;
00370 }
00371
00372 for( ion=0; ion <= nelem; ion++ )
00373 {
00374 for( j=0; j < (nelem + 2); j++ )
00375 {
00376 dense.xIonDense[nelem][j] = 0.;
00377 }
00378
00379 dense.xIonDense[nelem][ion] = 1.;
00380
00381 OpacityZero();
00382
00383
00384
00385 OpacityAdd1Element(nelem);
00386
00387
00388 strcpy( chFileName , elementnames.chElementNameShort[nelem] );
00389
00390
00391 strcat( chFileName , chNumbers[ion+1] );
00392
00393
00394 strcat( chFileName , ".opc" );
00395
00396
00397 ioFILENAME = open_data( chFileName, "w", AS_LOCAL_ONLY );
00398
00399 ilow = Heavy.ipHeavy[nelem][ion];
00400
00401 for( j=ilow-1; j < MIN2(rfield.nflux,continuum.KshellLimit); j++ )
00402 {
00403
00404 PrintE93( ioFILENAME , rfield.anu[j]*EVRYD );
00405 fprintf( ioFILENAME , "\t");
00406
00407
00408 PrintE93( ioFILENAME, (opac.opacity_abs[j]+opac.OpacStatic[j])/1e-18 );
00409 fprintf( ioFILENAME , "\n");
00410 }
00411
00412 fclose( ioFILENAME );
00413 }
00414
00415 prtPunOpacSummary();
00416 cdEXIT(EXIT_SUCCESS);
00417 }
00418
00419 return;
00420 }
00421
00422
00423 STATIC void prtPunOpacSummary(void)
00424 {
00425 fprintf(ioQQQ,"\n\nThe opacity files have been successfully created.\n");
00426 fprintf(ioQQQ,"The files have names that start with the first 4 characters of the element name.\n");
00427 fprintf(ioQQQ,"There is one file per ion and the number after the element name indicates the ion.\n");
00428 fprintf(ioQQQ,"The energies are in eV and the cross sections in megabarns.\n");
00429 fprintf(ioQQQ,"All end in \".opc\"\n");
00430 fprintf(ioQQQ,"The data only extend to the highest energy in this continuum source.\n");
00431 }