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 /*PrtContinuum print information about continuum if requested with PRINT CONTINUUM command, 00004 * called by PrtFinal */ 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 /* energies for the x-ray bands */ 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 /*print information about continuum if requested with PRINT CONTINUUM command */ 00040 /* this is number of ranges for adding x-ray bands */ 00041 00042 /* this stuff was always printed before version 84, and is now an option 00043 * to turn on information with PRINT CONTINUUM command */ 00044 00045 /* the default - not turned on, just return */ 00046 if( !prt.lgPrtCont ) 00047 { 00048 return; 00049 } 00050 00051 /* derive x-ray fluxes in various bands for later printout 00052 * 00053 * >>chng 97 jun 08, get rid of statement labels */ 00054 if( xraybd[0] < rfield.anu[rfield.nflux-1] ) 00055 { 00056 i1 = opac.ipCKshell - 10; 00057 /* get out to lower bound of first range */ 00058 while( (double)rfield.anu[i1-1] < xraybd[0] && i1 < rfield.nflux ) 00059 { 00060 i1 += 1; 00061 } 00062 /* now add up over that range */ 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[i1-1] + 00069 rfield.outlin[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 /* continuum not defined at x-ray energies, but zero it for safety */ 00078 for( i=0; i < NXBD; i++ ) 00079 { 00080 xbdsav[i] = 0.; 00081 } 00082 } 00083 00084 /* now print x-ray flux (photons s^-1, flux or lum) in various bands */ 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 /* renorm fLUX array before printing it */ 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[i] > 0. ) 00104 { 00105 fluxsv[i-iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]+1] = (realnum)((rfield.flux[i] +rfield.outlin[i] + rfield.outlin_noplot[i] + 00106 rfield.ConInterOut[i] )*radius.r1r0sq/ 00107 rfield.flux_total_incident[i]); 00108 } 00109 else 00110 { 00111 fluxsv[i-iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]+1] = 0.; 00112 } 00113 } 00114 /* renormalize whole continuum */ 00115 for( i=0; i < rfield.nflux; i++ ) 00116 { 00117 rfield.flux[i] = (realnum)(((rfield.flux[i] + 00118 rfield.ConInterOut[i]/rfield.anu[i])/rfield.widflx[i]+ rfield.outlin[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[i] + rfield.outlin[i] + rfield.outlin_noplot[i] +rfield.ConInterOut[i] ); 00145 } 00146 fprintf( ioQQQ, "\n" ); 00147 } 00148 return; 00149 }