00001 /* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*ParseCosmicRays parse the cosmic rays command */ 00004 #include "cddefines.h" 00005 #include "optimize.h" 00006 #include "hextra.h" 00007 #include "ionbal.h" 00008 #include "input.h" 00009 #include "parse.h" 00010 #include "parser.h" 00011 00012 /*ParseCosmicRays parse the cosmic rays command */ 00013 void ParseCosmicRays( Parser &p ) 00014 { 00015 int npar = 0; 00016 realnum a; 00017 double var; 00018 string ExtraPar; 00019 00020 DEBUG_ENTRY( "ParseCosmicRays()" ); 00021 00022 /* cosmic ray density, log of rate relative to background, log of H0 rate in neutral gas, 00023 * or density of rel. electrons, 00024 * quantity is log unless keyword linear appears */ 00025 /* if no number is present FFmtRead returns zero */ 00026 a = (realnum)p.FFmtRead(); 00027 if( p.lgEOL() ) 00028 a = 0.; 00029 00030 /* if keyword LINEAR not present, then log, and make linear */ 00031 if( !p.nMatch("LINE") ) 00032 a = (realnum)pow((realnum)10.f,a); 00033 /* a is now linear scale factor, or linear density, with default of 1 if no number */ 00034 00035 /* default is cosmic ray ionization rate relative to galactic background, but can 00036 * also give density, which was the only option originally */ 00037 if( p.nMatch("DENS") ) 00038 { 00039 if( p.lgEOL() ) 00040 { 00041 p.NoNumb("cosmic ray density"); 00042 } 00043 hextra.cryden = a; 00044 00045 /* optional power law density */ 00046 hextra.crpowr = (realnum)p.FFmtRead(); 00047 00048 /* option to specify a temp for non-rel electrons - but only when a density */ 00049 hextra.crtemp = (realnum)p.FFmtRead(); 00050 if( p.lgEOL() ) 00051 { 00052 /* relativistic limit (Balbus and McKee) */ 00053 hextra.crtemp = 2.6e9; 00054 } 00055 else 00056 { 00057 var = pow((realnum)10.f,hextra.crtemp); 00058 hextra.crtemp = (realnum)MIN2(var,2.6e9); 00059 } 00060 npar = 3; 00061 ExtraPar = "DENSITY"; 00062 } 00063 else if( p.nMatch( "RATE" ) ) 00064 { 00065 /* this sets rate - use stored density and rate for background to set 00066 * new density since code works with density */ 00067 ASSERT( a > 0. ); 00068 hextra.cryden = hextra.background_density * a / hextra.background_rate; 00069 hextra.crtemp = 2.6e9f; 00070 npar = 1; 00071 ExtraPar = "RATE"; 00072 } 00073 else if( p.nMatch( "BACKGROU" ) ) 00074 { 00075 /* >>chng 06 may 28, require explicit BACKGROUnd to hit background for safety */ 00076 /* cr relative to galactic background BACK - no check on string since default */ 00077 /* >>chng 04 mar 10, background is now 00078 * >>refer cr ion Williams, J.P., Bergin, E.A., Caseli, P., Myers, P.C., & Plume, R. 1998, ApJ, 503, 689 */ 00079 /* galactic background cosmic ray density to produce 00080 * secondary ionization rate quoted by Tielens and Hollenbach */ 00081 /* hextra.cryden = 2e-9f;*/ 00082 /* >>chng 99 jun 24, slight change to value 00083 * quoted by 00084 * >>refer cosmic ray ionization rate McKee, C.M., 1999, astro-ph 9901370 00085 * this will produce a total 00086 * secondary ionization rate of 2.5e-17 s^-1, as tested in 00087 * tsuite secondary.in. If each ionization produces 2.4 eV of heat, 00088 * the background heating rate should be 9.6e-29 * n*/ 00089 /* >>chng 00 nov 28, changed density to 4.9e-9 to reproduce TH85a 00090 * when photoionization is turned off. 00091 >>refer cosmic ray ionization rate Tielens, A.G.G.M., & Hollenbach, D., 1998, ApJ, 291, 722 00092 */ 00093 /* hextra.cryden = 7.07e-9f;*/ 00094 /* this value reproduces the TH cr ionization rate when the factor 00095 * of 0.46 is included. This will directly go onto the h ionization rate 00096 * without the factor of 0.46 there. this is necessary for the more 00097 * general case where cr ionization is actually self-consistently determined 00098 * from rate hot electrons injected into the plasma */ 00099 /*hextra.cryden = 2.25e-9f;*/ 00100 ASSERT( a > 0. ); 00101 hextra.cryden = hextra.background_density * a; 00102 hextra.crtemp = 2.6e9f; 00103 npar = 1; 00104 ExtraPar = "BACKGROUND"; 00105 } 00106 else if( p.nMatch( "EQUI" ) ) 00107 { 00108 /* equipartition cosmic rays, set from B */ 00109 hextra.lg_CR_B_equipartition = true; 00110 /* this has to be positive for cr's to be on 00111 * it will be reevaluated when B is known */ 00112 hextra.cryden = SMALLFLOAT; 00113 hextra.crtemp = 2.6e9f; 00114 } 00115 00116 else 00117 { 00118 /* no keyword found */ 00119 fprintf( ioQQQ, " There must be a keyword on this COSMIC RAY command.\n" ); 00120 fprintf( ioQQQ, " The keywords are DENSITY, RATE, BACKGROUND, and EQUIPARTITION.\n" ); 00121 cdEXIT(EXIT_FAILURE); 00122 } 00123 00124 /* this is current cosmic ray density divided by background - used in 00125 * a few chemical reactions */ 00126 hextra.cryden_ov_background = hextra.cryden / hextra.background_density; 00127 /* >>chng 05 jan 05, 00128 * set the cr ionization rate to very rough value, before we have enough 00129 * information to evaluate it - may be needed in initial guess of H and He ionization*/ 00130 ionbal.CosRayIonRate = hextra.cryden_ov_background * 2.5e-17; 00131 00132 /* vary option */ 00133 if( optimize.lgVarOn && ExtraPar.length() > 0 ) 00134 { 00135 /* will be one parameter */ 00136 optimize.nvarxt[optimize.nparm] = npar; 00137 sprintf( optimize.chVarFmt[optimize.nparm], "COSMic rays %s= %%f LOG", ExtraPar.c_str() ); 00138 /* log of cosmic rays rates relative to background */ 00139 optimize.vparm[0][optimize.nparm] = (realnum)log10(a); 00140 if( npar == 3 ) 00141 { 00142 strcat( optimize.chVarFmt[optimize.nparm], " %f %f" ); 00143 optimize.vparm[1][optimize.nparm] = hextra.crpowr; 00144 optimize.vparm[2][optimize.nparm] = realnum(log10(hextra.crtemp)); 00145 } 00146 /* array index for where to write */ 00147 optimize.nvfpnt[optimize.nparm] = input.nRead; 00148 /* the increment in the first steps away from the original value */ 00149 optimize.vincr[optimize.nparm] = 0.2f; 00150 ++optimize.nparm; 00151 } 00152 00153 return; 00154 }