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 /*ParseConstant parse parameters from the 'constant ...' command */ 00004 #include "cddefines.h" 00005 #include "physconst.h" 00006 #include "optimize.h" 00007 #include "thermal.h" 00008 #include "dense.h" 00009 #include "pressure.h" 00010 #include "phycon.h" 00011 #include "input.h" 00012 #include "parse.h" 00013 00014 void ParseConstant(char *chCard ) 00015 { 00016 bool lgEOL; 00017 long int i; 00018 00019 DEBUG_ENTRY( "ParseConstant()" ); 00020 00021 if( nMatch("GRAI",chCard) && nMatch("TEMP",chCard) ) 00022 { 00023 /* constant grain temperature command */ 00024 i = 5; 00025 thermal.ConstGrainTemp = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00026 00027 /* if linear option is not on the line, convert to exponent if <= 10 */ 00028 if( !nMatch("LINE",chCard) ) 00029 { 00030 if( thermal.ConstGrainTemp <= 10. ) 00031 thermal.ConstGrainTemp = (realnum)pow((realnum)10.f,thermal.ConstGrainTemp); 00032 } 00033 00034 if( lgEOL ) 00035 { 00036 NoNumb(chCard); 00037 } 00038 } 00039 00040 else if( nMatch("TEMP",chCard) ) 00041 { 00042 /* a constant temperature model */ 00043 thermal.lgTemperatureConstant = true; 00044 thermal.lgTemperatureConstantCommandParsed = true; 00045 00046 /* this is an option to specify the temperature in different units 00047 * keV, eV now supported */ 00048 realnum convert_to_Kelvin = 1; 00049 if( nMatch(" EV ",chCard) ) 00050 { 00051 convert_to_Kelvin = (realnum)EVDEGK; 00052 } 00053 else if( nMatch(" KEV",chCard) ) 00054 { 00055 convert_to_Kelvin = (realnum)(EVDEGK * 1000.); 00056 } 00057 00058 i = 5; 00059 /* this is the "force" temperature. same var used in force temp 00060 * command, but lgTSetOn is not set so then allowed to vary 00061 * so constant temperature requires both lgTSetOn true and ConstTemp > 0 */ 00062 thermal.ConstTemp = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00063 if( lgEOL ) 00064 NoNumb(chCard); 00065 00066 /* if linear option is not on the line and T<=10, assume number is log */ 00067 if( !nMatch("LINE",chCard) ) 00068 { 00069 if( thermal.ConstTemp <= 10. ) 00070 thermal.ConstTemp = (realnum)pow((realnum)10.f,thermal.ConstTemp); 00071 } 00072 /* do units conversion here */ 00073 thermal.ConstTemp *= convert_to_Kelvin; 00074 00075 /* check that temperature is not less than 3K */ 00076 if( thermal.ConstTemp < 3. ) 00077 { 00078 fprintf( ioQQQ, " PROBLEM TE<3K, reset to 3K.\n" ); 00079 thermal.ConstTemp = 3.; 00080 } 00081 00082 /* bail if temperature is too high */ 00083 /* >>chng 06 feb 14, can cannot do a temperature of 1e10 - gaunt factors overrun 00084 * energy array. temp must be less than 1e10 00085 * >>chng 06 feb 27, change 1e10 to phycon.TEMP_LIMIT_HIGH - which is the same */ 00086 if( thermal.ConstTemp > phycon.TEMP_LIMIT_HIGH ) 00087 { 00088 fprintf( ioQQQ, " The constant temperature (%.3e K) is > %.3e K, Cloudy cannot handle this temperature.\n", 00089 thermal.ConstTemp, 00090 phycon.TEMP_LIMIT_HIGH); 00091 fprintf( ioQQQ, " The upper temperature limit is %.3e K.\n" , 00092 phycon.TEMP_LIMIT_HIGH); 00093 fprintf( ioQQQ, " Sorry.\n" ); 00094 cdEXIT(EXIT_FAILURE); 00095 } 00096 00097 /* set the real electron temperature to the forced value */ 00098 TempChange(thermal.ConstTemp , false); 00099 00100 /* vary option */ 00101 if( optimize.lgVarOn ) 00102 { 00103 /* no luminosity options on vary */ 00104 optimize.nvarxt[optimize.nparm] = 1; 00105 strcpy( optimize.chVarFmt[optimize.nparm], "CONStant TEMP %f" ); 00106 00107 /* pointer to where to write */ 00108 optimize.nvfpnt[optimize.nparm] = input.nRead; 00109 00110 /* log of temp will be pointer */ 00111 optimize.vparm[0][optimize.nparm] = (realnum)log10(phycon.te); 00112 optimize.vincr[optimize.nparm] = 0.1f; 00113 ++optimize.nparm; 00114 } 00115 } 00116 00117 else if( nMatch("DENS",chCard) ) 00118 { 00119 /* constant density */ 00120 strcpy( dense.chDenseLaw, "CDEN" ); 00121 /* turn off radiation pressure */ 00122 pressure.lgPres_radiation_ON = false; 00123 pressure.lgPres_magnetic_ON = false; 00124 pressure.lgPres_ram_ON = false; 00125 } 00126 00127 else if( nMatch("PRES",chCard) ) 00128 { 00129 /* constant pressure */ 00130 strcpy( dense.chDenseLaw, "CPRE" ); 00131 00132 /* >>chng 06 jun 20, add reset option, to reset density to keep 00133 * initial pressure itself constant from iteration to iteration, 00134 * rather than initial density */ 00135 if( nMatch("RESE",chCard) ) 00136 { 00137 /* this says not to keep initial density constant, 00138 * reset it to keep pressure const */ 00139 dense.lgDenseInitConstant = false; 00140 } 00141 else 00142 { 00143 /* this is default, says keep initial density constant, 00144 * so pressure from iter to iter not really const */ 00145 dense.lgDenseInitConstant = true; 00146 } 00147 00148 if( nMatch(" GAS",chCard) ) 00149 { 00150 /* constant gas pressure (no radiation) 00151 * turn off radiation pressure */ 00152 pressure.lgPres_radiation_ON = false; 00153 00154 /* turn off incident continuum */ 00155 pressure.lgContRadPresOn = false; 00156 00157 /* turn off magnetic and ram pressure */ 00158 pressure.lgPres_magnetic_ON = false; 00159 pressure.lgPres_ram_ON = false; 00160 00161 /* optional second number is power law index */ 00162 i = 4; 00163 pressure.PresPowerlaw = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00164 } 00165 00166 else 00167 { 00168 /* constant total pressure, gas+rad+incident continuum 00169 * turn on radiation pressure */ 00170 pressure.lgPres_radiation_ON = true; 00171 pressure.lgPres_magnetic_ON = true; 00172 pressure.lgPres_ram_ON = true; 00173 00174 /* option to turn off continuum pressure */ 00175 if( nMatch("O CO",chCard) ) 00176 { 00177 pressure.lgContRadPresOn = false; 00178 } 00179 else 00180 { 00181 pressure.lgContRadPresOn = true; 00182 } 00183 00184 /* option to not abort when too much radiation pressure, no abort */ 00185 if( nMatch("O AB",chCard) ) 00186 { 00187 pressure.lgRadPresAbortOK = false; 00188 } 00189 else 00190 { 00191 pressure.lgRadPresAbortOK = true; 00192 } 00193 /* there is no optional power law option for constant total pressure */ 00194 pressure.PresPowerlaw = 0.; 00195 00196 /* option to set pressure */ 00197 if( nMatch(" SET",chCard) ) 00198 { 00199 /* number on line is log of nT - option to specify initial pressure */ 00200 pressure.lgPressureInitialSpecified = true; 00201 /* this is log of nT product - if not present then set zero */ 00202 i = 5; 00203 pressure.PressureInitialSpecified = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00204 if( lgEOL ) 00205 NoNumb( chCard ); 00206 else 00207 /* pressure in nkT units */ 00208 pressure.PressureInitialSpecified = pow(10.,pressure.PressureInitialSpecified ) * 00209 BOLTZMANN; 00210 } 00211 else 00212 pressure.lgPressureInitialSpecified = false; 00213 } 00214 00215 /* vary option */ 00216 if( optimize.lgVarOn &&pressure.lgPressureInitialSpecified ) 00217 { 00218 /* no options on vary */ 00219 optimize.nvarxt[optimize.nparm] = 1; 00220 strcpy( optimize.chVarFmt[optimize.nparm], "CONStant PRESsure SET %f" ); 00221 00222 /* pointer to where to write */ 00223 optimize.nvfpnt[optimize.nparm] = input.nRead; 00224 00225 /* log of temp will be pointer */ 00226 optimize.vparm[0][optimize.nparm] = (realnum)log10(pressure.PressureInitialSpecified); 00227 optimize.vincr[optimize.nparm] = 0.1f; 00228 ++optimize.nparm; 00229 } 00230 } 00231 00232 else 00233 { 00234 /* no keys were recognized */ 00235 fprintf( ioQQQ, " The keyword should be TEMPerature, DENSity, GAS or PRESsure, sorry.\n" ); 00236 cdEXIT(EXIT_FAILURE); 00237 } 00238 return; 00239 }