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 realnum 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 return;
00049 }
00050
00051
00052
00053
00054 if( xraybd[0] < rfield.anu[rfield.nflux-1] )
00055 {
00056 i1 = opac.ipCKshell - 10;
00057
00058 while( (double)rfield.anu[i1-1] < xraybd[0] && i1 < rfield.nflux )
00059 {
00060 i1 += 1;
00061 }
00062
00063 for( i=0; i < NXBD; i++ )
00064 {
00065 xbdsav[i] = 0.;
00066 while( (double)rfield.anu[i1-1] < xraybd[i+1] && i1 < rfield.nflux )
00067 {
00068 xbdsav[i] += rfield.flux[0][i1-1] +
00069 rfield.outlin[0][i1-1] +rfield.outlin_noplot[i1-1] + rfield.ConInterOut[i1-1];
00070 i1 += 1;
00071 }
00072 xbdsav[i] *= (realnum)radius.r1r0sq;
00073 }
00074 }
00075 else
00076 {
00077
00078 for( i=0; i < NXBD; i++ )
00079 {
00080 xbdsav[i] = 0.;
00081 }
00082 }
00083
00084
00085 if( xbdsav[0] > 0 )
00086 {
00087 fprintf( ioQQQ, "\n 0.1-0.5kev:" );
00088 for(i=0; i < NXBD; i++)
00089 fprintf( ioQQQ, "%8.2e", xbdsav[i] );
00090 fprintf( ioQQQ, " 0.5-1.0kev:\n" );
00091 }
00092
00093
00094 if( iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2] + 1 > NFLUXSV )
00095 {
00096 fprintf( ioQQQ, " PCONTN: not enough cells in flux_total_incident, need%4ld\n",
00097 iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s] - iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2] + 1 );
00098 cdEXIT(EXIT_FAILURE);
00099 }
00100
00101 for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1; i < iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]; i++ )
00102 {
00103 if( rfield.flux_total_incident[0][i] > 0. )
00104 {
00105 fluxsv[i-iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]+1] = (realnum)((rfield.flux[0][i] +rfield.outlin[0][i] + rfield.outlin_noplot[i] +
00106 rfield.ConInterOut[i] )*radius.r1r0sq/
00107 rfield.flux_total_incident[0][i]);
00108 }
00109 else
00110 {
00111 fluxsv[i-iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]+1] = 0.;
00112 }
00113 }
00114
00115 for( i=0; i < rfield.nflux; i++ )
00116 {
00117 rfield.flux[0][i] = (realnum)(((rfield.flux[0][i] +
00118 rfield.ConInterOut[i]/rfield.anu[i])/rfield.widflx[i]+ rfield.outlin[0][i] + rfield.outlin_noplot[i])*
00119 radius.r1r0sq);
00120 }
00121
00122 fprintf( ioQQQ,
00123 "\n\n Normalised continuum\n" );
00124
00125 for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]; i <= iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]; i += 3 )
00126 {
00127 fprintf( ioQQQ, "%7.3f%6.3f", rfield.anu[i-1], fluxsv[i-iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]] );
00128 }
00129 fprintf( ioQQQ, "\n" );
00130
00131 fprintf( ioQQQ,
00132 "\n Emergent continuum - phot/ryd/cm2/s (r in)\n" );
00133
00134 nline = ((rfield.nflux - iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2] - 1)/4 + 7)/7;
00135 ninc = nline*4;
00136 for( j=0; j < nline; j++ )
00137 {
00138 i1 = j*4 + iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2];
00139 fprintf( ioQQQ, " " );
00140
00141 for( i=i1; i<rfield.nflux; i = i + ninc)
00142 {
00143 fprintf( ioQQQ, "%6.3f%10.2e", rfield.anu[i],
00144 rfield.flux[0][i] + rfield.outlin[0][i] + rfield.outlin_noplot[i] +rfield.ConInterOut[i] );
00145 }
00146 fprintf( ioQQQ, "\n" );
00147 }
00148 return;
00149 }