00001
00002
00003
00004
00005 #include "cddefines.h"
00006 #include "rfield.h"
00007 #include "iso.h"
00008 #include "radius.h"
00009 #include "opacity.h"
00010 #include "prt.h"
00011
00012 static const int NFLUXSV = 360;
00013 static const int NXBD = 9;
00014
00015 void PrtContinuum(void)
00016 {
00017 long int i,
00018 i1,
00019 j,
00020 ninc,
00021 nline;
00022 float fluxsv[NFLUXSV],
00023 xbdsav[NXBD];
00024
00025
00026 static double xraybd[NXBD + 1]={
00027 7.3498,
00028 36.8,
00029 73.5,
00030 110.3,
00031 147.1,
00032 183.8,
00033 220.6,
00034 367.6,
00035 551.5,
00036 735.3};
00037
00038 DEBUG_ENTRY( "PrtContinuum()" );
00039
00040
00041
00042
00043
00044
00045
00046 if( !prt.lgPrtCont )
00047 {
00048 DEBUG_EXIT( "PrtContinuum()" );
00049 return;
00050 }
00051
00052
00053
00054
00055 if( xraybd[0] < rfield.anu[rfield.nflux-1] )
00056 {
00057 i1 = opac.ipCKshell - 10;
00058
00059 while( (double)rfield.anu[i1-1] < xraybd[0] && i1 < rfield.nflux )
00060 {
00061 i1 += 1;
00062 }
00063
00064 for( i=0; i < NXBD; i++ )
00065 {
00066 xbdsav[i] = 0.;
00067 while( (double)rfield.anu[i1-1] < xraybd[i+1] && i1 < rfield.nflux )
00068 {
00069 xbdsav[i] += rfield.flux[i1-1] +
00070 rfield.outlin[i1-1] +rfield.outlin_noplot[i1-1] + rfield.ConInterOut[i1-1];
00071 i1 += 1;
00072 }
00073 xbdsav[i] *= (float)radius.r1r0sq;
00074 }
00075 }
00076 else
00077 {
00078
00079 for( i=0; i < NXBD; i++ )
00080 {
00081 xbdsav[i] = 0.;
00082 }
00083 }
00084
00085
00086 if( xbdsav[0] > 0 )
00087 {
00088 fprintf( ioQQQ, "\n 0.1-0.5kev:" );
00089 for(i=0; i < NXBD; i++)
00090 fprintf( ioQQQ, "%8.2e", xbdsav[i] );
00091 fprintf( ioQQQ, " 0.5-1.0kev:\n" );
00092 }
00093
00094
00095 if( iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2] + 1 > NFLUXSV )
00096 {
00097 fprintf( ioQQQ, " PCONTN: not enough cells in FluxSave, need%4ld\n",
00098 iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2] + 1 );
00099 puts( "[Stop in PrtContinuum]" );
00100 cdEXIT(EXIT_FAILURE);
00101 }
00102
00103 for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1; i < iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]; i++ )
00104 {
00105 if( rfield.FluxSave[i] > 0. )
00106 {
00107 fluxsv[i-iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]+1] = (float)((rfield.flux[i] +rfield.outlin[i] + rfield.outlin_noplot[i] +
00108 rfield.ConInterOut[i] )*radius.r1r0sq/
00109 rfield.FluxSave[i]);
00110 }
00111 else
00112 {
00113 fluxsv[i-iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]+1] = 0.;
00114 }
00115 }
00116
00117 for( i=0; i < rfield.nflux; i++ )
00118 {
00119 rfield.flux[i] = (float)(((rfield.flux[i] +
00120 rfield.ConInterOut[i]/rfield.anu[i])/rfield.widflx[i]+ rfield.outlin[i] + rfield.outlin_noplot[i])*
00121 radius.r1r0sq);
00122 }
00123
00124 fprintf( ioQQQ,
00125 "\n\n Normalised continuum\n" );
00126
00127 for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]; i <= iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]; i += 3 )
00128 {
00129 fprintf( ioQQQ, "%7.3f%6.3f", rfield.anu[i-1], fluxsv[i-iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]] );
00130 }
00131 fprintf( ioQQQ, "\n" );
00132
00133 fprintf( ioQQQ,
00134 "\n Emergent continuum - phot/ryd/cm2/s (r in)\n" );
00135
00136 nline = ((rfield.nflux - iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2] - 1)/4 + 7)/7;
00137 ninc = nline*4;
00138 for( j=0; j < nline; j++ )
00139 {
00140 i1 = j*4 + iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2];
00141 fprintf( ioQQQ, " " );
00142
00143 for( i=i1; i<rfield.nflux; i = i + ninc)
00144 {
00145 fprintf( ioQQQ, "%6.3f%10.2e", rfield.anu[i],
00146 rfield.flux[i] + rfield.outlin[i] + rfield.outlin_noplot[i] +rfield.ConInterOut[i] );
00147 }
00148 fprintf( ioQQQ, "\n" );
00149 }
00150
00151 DEBUG_EXIT( "PrtContinuum()" );
00152 return;
00153 }