00001
00002
00003
00004 #include "cddefines.h"
00005 #include "energy.h"
00006 #include "continuum.h"
00007 #include "rfield.h"
00008 #include "doppvel.h"
00009 #include "geometry.h"
00010 #include "hextra.h"
00011 #include "ipoint.h"
00012 #include "lines.h"
00013 #include "lines_service.h"
00014 #include "opacity.h"
00015 #include "physconst.h"
00016 #include "prt.h"
00017 #include "radius.h"
00018 #include "secondaries.h"
00019 #include "struc.h"
00020 #include "thermal.h"
00021 #include "warnings.h"
00022 #include "wind.h"
00023 #include "dynamics.h"
00024
00025 static const char *ENERGY_RYD = "Ryd";
00026 static const char *ENERGY_ERG = "erg";
00027 static const char *ENERGY_EV = "eV";
00028 static const char *ENERGY_KEV = "keV";
00029 static const char *ENERGY_MEV = "MeV";
00030 static const char *ENERGY_WN = "cm^-1";
00031 static const char *ENERGY_CM = "cm";
00032 static const char *ENERGY_MM = "mm";
00033 static const char *ENERGY_MICRON = "um";
00034 static const char *ENERGY_NM = "nm";
00035 static const char *ENERGY_A = "A";
00036 static const char *ENERGY_HZ = "Hz";
00037 static const char *ENERGY_KHZ = "kHz";
00038 static const char *ENERGY_MHZ = "MHz";
00039 static const char *ENERGY_GHZ = "GHz";
00040 static const char *ENERGY_K = "K";
00041
00042 inline bool isSameUnit(const char *unit, const char *ref)
00043 {
00044 return strcmp(unit,ref) == 0;
00045 }
00046
00047 const char *StandardEnergyUnit(const char *chCard)
00048 {
00049 DEBUG_ENTRY( "StandardEnergyUnit()" );
00050
00051
00052 if( nMatch(" MIC",chCard) )
00053 {
00054
00055 return ENERGY_MICRON;
00056 }
00057 else if( nMatch(" EV ",chCard) )
00058 {
00059
00060 return ENERGY_EV;
00061 }
00062 else if( nMatch(" KEV",chCard) )
00063 {
00064
00065 return ENERGY_KEV;
00066 }
00067 else if( nMatch(" MEV",chCard) )
00068 {
00069
00070 return ENERGY_MEV;
00071 }
00072 else if( nMatch("WAVE",chCard) )
00073 {
00074
00075 return ENERGY_WN;
00076 }
00077 else if( nMatch("CENT",chCard) || nMatch(" CM ",chCard) )
00078 {
00079
00080 return ENERGY_CM;
00081 }
00082 else if( nMatch(" MM ",chCard) )
00083 {
00084
00085 return ENERGY_MM;
00086 }
00087 else if( nMatch(" NM ",chCard) )
00088 {
00089
00090 return ENERGY_NM;
00091 }
00092 else if( nMatch("ANGS",chCard) )
00093 {
00094
00095 return ENERGY_A;
00096 }
00097 else if( nMatch(" HZ ",chCard) )
00098 {
00099
00100 return ENERGY_HZ;
00101 }
00102 else if( nMatch(" KHZ",chCard) )
00103 {
00104
00105 return ENERGY_KHZ;
00106 }
00107 else if( nMatch(" MHZ",chCard) )
00108 {
00109
00110 return ENERGY_MHZ;
00111 }
00112 else if( nMatch(" GHZ",chCard) )
00113 {
00114
00115 return ENERGY_GHZ;
00116 }
00117 else if( nMatch("KELV",chCard) || nMatch(" K ",chCard) )
00118 {
00119
00120 return ENERGY_K;
00121 }
00122
00123 else if( nMatch("ERG ",chCard) )
00124 {
00125
00126 return ENERGY_ERG;
00127 }
00128 else if( nMatch(" RYD",chCard) )
00129 {
00130
00131 return ENERGY_RYD;
00132 }
00133 else
00134 {
00135 fprintf( ioQQQ, " No energy / wavelength unit was recognized on this line:\n %s\n\n", chCard );
00136 fprintf( ioQQQ, " See Hazy for details.\n" );
00137 cdEXIT(EXIT_FAILURE);
00138 }
00139 }
00140
00141 double Energy::get(const char *unit) const
00142 {
00143 DEBUG_ENTRY( "Energy::get()" );
00144
00145
00146 if( isSameUnit(unit,ENERGY_RYD) )
00147 {
00148 return Ryd();
00149 }
00150 else if( isSameUnit(unit,ENERGY_ERG) )
00151 {
00152 return Erg();
00153 }
00154 else if( isSameUnit(unit,ENERGY_EV) )
00155 {
00156 return eV();
00157 }
00158 else if( isSameUnit(unit,ENERGY_KEV) )
00159 {
00160 return keV();
00161 }
00162 else if( isSameUnit(unit,ENERGY_MEV) )
00163 {
00164 return MeV();
00165 }
00166 else if( isSameUnit(unit,ENERGY_WN) )
00167 {
00168 return WN();
00169 }
00170 else if( isSameUnit(unit,ENERGY_CM) )
00171 {
00172 return cm();
00173 }
00174 else if( isSameUnit(unit,ENERGY_MM) )
00175 {
00176 return mm();
00177 }
00178 else if( isSameUnit(unit,ENERGY_MICRON) )
00179 {
00180 return micron();
00181 }
00182 else if( isSameUnit(unit,ENERGY_NM) )
00183 {
00184 return nm();
00185 }
00186 else if( isSameUnit(unit,ENERGY_A) )
00187 {
00188 return Angstrom();
00189 }
00190 else if( isSameUnit(unit,ENERGY_HZ) )
00191 {
00192 return Hz();
00193 }
00194 else if( isSameUnit(unit,ENERGY_KHZ) )
00195 {
00196 return kHz();
00197 }
00198 else if( isSameUnit(unit,ENERGY_MHZ) )
00199 {
00200 return MHz();
00201 }
00202 else if( isSameUnit(unit,ENERGY_GHZ) )
00203 {
00204 return GHz();
00205 }
00206 else if( isSameUnit(unit,ENERGY_K) )
00207 {
00208 return K();
00209 }
00210 else
00211 {
00212 fprintf( ioQQQ, " insane units in Energy::get: \"%s\"\n", unit );
00213 cdEXIT(EXIT_FAILURE);
00214 }
00215 }
00216
00217 void Energy::set(double value, const char *unit)
00218 {
00219 DEBUG_ENTRY( "Energy::set()" );
00220
00221
00222 if( isSameUnit(unit,ENERGY_RYD) )
00223 {
00224 set(value);
00225 }
00226 else if( isSameUnit(unit,ENERGY_ERG) )
00227 {
00228 set(value/EN1RYD);
00229 }
00230 else if( isSameUnit(unit,ENERGY_MEV) )
00231 {
00232 set(1e6*value/EVRYD);
00233 }
00234 else if( isSameUnit(unit,ENERGY_KEV) )
00235 {
00236 set(1e3*value/EVRYD);
00237 }
00238 else if( isSameUnit(unit,ENERGY_EV) )
00239 {
00240 set(value/EVRYD);
00241 }
00242 else if( isSameUnit(unit,ENERGY_WN) )
00243 {
00244 set(value/RYD_INF);
00245 }
00246 else if( isSameUnit(unit,ENERGY_A) )
00247 {
00248 set(RYDLAM/value);
00249 }
00250 else if( isSameUnit(unit,ENERGY_NM) )
00251 {
00252 set(RYDLAM/(1.e1*value));
00253 }
00254 else if( isSameUnit(unit,ENERGY_MICRON) )
00255 {
00256 set(RYDLAM/(1.e4*value));
00257 }
00258 else if( isSameUnit(unit,ENERGY_MM) )
00259 {
00260 set(RYDLAM/(1.e7*value));
00261 }
00262 else if( isSameUnit(unit,ENERGY_CM) )
00263 {
00264 set(RYDLAM/(1.e8*value));
00265 }
00266 else if( isSameUnit(unit,ENERGY_HZ) )
00267 {
00268 set(value/FR1RYD);
00269 }
00270 else if( isSameUnit(unit,ENERGY_KHZ) )
00271 {
00272 set(1e3*value/FR1RYD);
00273 }
00274 else if( isSameUnit(unit,ENERGY_MHZ) )
00275 {
00276 set(1e6*value/FR1RYD);
00277 }
00278 else if( isSameUnit(unit,ENERGY_GHZ) )
00279 {
00280 set(1e9*value/FR1RYD);
00281 }
00282 else if( isSameUnit(unit,ENERGY_K) )
00283 {
00284 set(value/TE1RYD);
00285 }
00286 else
00287 {
00288 fprintf( ioQQQ, " insane units in Energy::set: \"%s\"\n", unit );
00289 cdEXIT(EXIT_FAILURE);
00290 }
00291 }
00292
00293 void EnergyEntry::p_set_ip()
00294 {
00295 DEBUG_ENTRY( "EnergyEntry::p_set_ip()" );
00296
00297 double energy = Ryd();
00298 if( energy < rfield.emm || energy > rfield.egamry )
00299 {
00300 fprintf( ioQQQ, " The energy %g Ryd is not within the default Cloudy range\n", energy );
00301 cdEXIT(EXIT_FAILURE);
00302 }
00303 p_ip = ipoint(energy) - 1;
00304 ASSERT( p_ip >= 0 );
00305 }
00306
00307 STATIC double sum_radiation( void );
00308
00309
00310 STATIC void badprt(
00311
00312 double total);
00313
00314 bool lgConserveEnergy()
00315 {
00316 double outin,
00317 outout,
00318 outtot,
00319 sum_coolants,
00320 sum_recom_lines,
00321 flow_energy;
00322
00323 char chLine[INPUT_LINE_LENGTH];
00324
00325 DEBUG_ENTRY( "lgConserveEnergy()" );
00326
00327
00328
00329 outsum(&outtot,&outin,&outout);
00330
00331 sum_recom_lines = totlin('r');
00332 if( sum_recom_lines == 0. )
00333 {
00334 sprintf( chLine, " !Total recombination line energy is 0." );
00335 bangin(chLine);
00336 }
00337
00338
00339 sum_coolants = totlin('c');
00340 if( sum_coolants == 0. )
00341 {
00342 sprintf( chLine, " !Total cooling is zero." );
00343 bangin(chLine);
00344 }
00345
00346 if( !(wind.lgBallistic() || wind.lgStatic()) )
00347 {
00348
00349
00350 flow_energy = (2.5*struc.GasPressure[0]+0.5*struc.DenMass[0]*wind.windv0*wind.windv0)
00351 *(-wind.windv0);
00352 }
00353 else
00354 {
00355 flow_energy = 0.;
00356 }
00357
00358 double sum_rad = sum_radiation();
00359
00360 if( 0 && geometry.lgSphere && (!thermal.lgTemperatureConstant) && (!dynamics.lgTimeDependentStatic) &&
00361 (hextra.cryden == 0.) && ((hextra.TurbHeat+DoppVel.DispScale) == 0.) && !secondaries.lgCSetOn )
00362 {
00363 if( fabs( 1. - sum_rad/(continuum.TotalLumin*POW2(radius.rinner)) ) > 0.01 )
00364 fprintf( ioQQQ, "PROBLEM DISASTER lgConserveEnergy failed %li\t%li\t%e\t%e\t%e\n",
00365 iteration, nzone, sum_rad, continuum.TotalLumin*POW2(radius.rinner),
00366 sum_rad/(continuum.TotalLumin*POW2(radius.rinner)) );
00367 }
00368
00369
00370
00371
00372
00373
00374 if( sum_coolants + sum_recom_lines + flow_energy > continuum.TotalLumin*geometry.covgeo &&
00375 !thermal.lgTemperatureConstant && geometry.iEmissPower == 2 )
00376 {
00377 if( (hextra.cryden == 0.) && ((hextra.TurbHeat+DoppVel.DispScale) == 0.) && !secondaries.lgCSetOn )
00378 {
00379
00380 sprintf( chLine,
00381 " W-Radiated luminosity (cool+rec+wind=%.2e+%.2e+%.2e) is greater than that in incident radiation (TotalLumin=%8.2e). Power radiated was %8.2e",
00382 sum_coolants, sum_recom_lines, flow_energy , continuum.TotalLumin, thermal.power );
00383 warnin(chLine);
00384
00385 fprintf( ioQQQ, "\n\n DISASTER This calculation DID NOT CONSERVE ENERGY!\n\n\n" );
00386
00387 lgAbort = true;
00388
00389
00390 if( opac.lgCaseB )
00391 fprintf( ioQQQ, "\n The CASE B command was entered - this can have unphysical effects, including non-conservation of energy.\n Why was it needed?\n\n" );
00392
00393
00394 badprt(continuum.TotalLumin);
00395
00396 sprintf( chLine, " W-Something is really wrong: the ratio of radiated to incident luminosity is %.2e.",
00397 (sum_coolants + sum_recom_lines)/continuum.TotalLumin );
00398 warnin(chLine);
00399
00400
00401 if( thermal.ConstTemp > 0. )
00402 {
00403 fprintf( ioQQQ, " This may have been caused by the FORCE TE command.\n" );
00404 fprintf( ioQQQ, " Remove it and run again.\n" );
00405 }
00406 else
00407 return false;
00408 }
00409 }
00410
00411 return true;
00412 }
00413
00414 STATIC double sum_radiation( void )
00415 {
00416 double flxref = 0., flxatt = 0., conem = 0.;
00417 long nEmType = 0;
00418 double sum = 0.;
00419
00420 for( long j=0; j<rfield.nflux; j++ )
00421 {
00422
00423 flxref += (rfield.ConRefIncid[nEmType][j] + rfield.ConEmitReflec[nEmType][j] + rfield.reflin[nEmType][j]) * rfield.anu[j];
00424 ASSERT( flxref >= 0. );
00425
00426
00427 flxatt += rfield.flux[nEmType][j]*radius.r1r0sq*rfield.trans_coef_total[j] * rfield.anu[j];
00428 ASSERT( flxatt >= 0. );
00429
00430
00431 conem += (rfield.ConEmitOut[nEmType][j] + rfield.outlin[nEmType][j] + rfield.outlin_noplot[j])*radius.r1r0sq*geometry.covgeo * rfield.anu[j];
00432 ASSERT( conem >= 0. );
00433 }
00434
00435 sum = (flxref + flxatt + conem) * EN1RYD * POW2(radius.rinner);
00436 ASSERT( sum >= 0. );
00437
00438 return sum;
00439 }
00440
00441
00442 STATIC void badprt(
00443
00444 double total)
00445 {
00446
00447 static const double RATIO = 0.02;
00448 char chInfo;
00449 long int i;
00450 realnum sum_coolants,
00451 sum_recom_lines;
00452
00453 DEBUG_ENTRY( "badprt()" );
00454
00455 fprintf( ioQQQ, " badprt: all entries with greater than%6.2f%% of incident continuum follow.\n",
00456 RATIO*100. );
00457 fprintf( ioQQQ, " badprt: Intensities are relative to total energy in incident continuum.\n" );
00458
00459
00460 chInfo = 'r';
00461 sum_recom_lines = (realnum)totlin('r');
00462 fprintf( ioQQQ,
00463 " Sum of energy in recombination lines is %.2e relative to total incident radiation is %.2e\n",
00464 sum_recom_lines,
00465 sum_recom_lines/MAX2(1e-30,total) );
00466
00467 fprintf(ioQQQ," all strong information lines \n line wl ener/total\n");
00468
00469 for( i=0; i < LineSave.nsum; i++ )
00470 {
00471 if( LineSv[i].chSumTyp == chInfo && LineSv[i].SumLine[0]/total > RATIO )
00472 {
00473 fprintf( ioQQQ, " %4.4s ", LineSv[i].chALab );
00474 prt_wl( ioQQQ, LineSv[i].wavelength );
00475 fprintf( ioQQQ, " %7.3f %c\n", LineSv[i].SumLine[0]/total, chInfo );
00476 }
00477 }
00478
00479 fprintf(ioQQQ," all strong cooling lines \n line wl ener/total\n");
00480 chInfo = 'c';
00481 sum_coolants = (realnum)totlin('c');
00482 fprintf( ioQQQ, " Sum of coolants (abs) = %.2e (rel)= %.2e\n",
00483 sum_coolants, sum_coolants/MAX2(1e-30,total) );
00484 for( i=0; i < LineSave.nsum; i++ )
00485 {
00486 if( LineSv[i].chSumTyp == chInfo && LineSv[i].SumLine[0]/total > RATIO )
00487 {
00488 fprintf( ioQQQ, " %4.4s ", LineSv[i].chALab );
00489 prt_wl( ioQQQ, LineSv[i].wavelength );
00490 fprintf( ioQQQ, " %7.3f %c\n", LineSv[i].SumLine[0]/total, chInfo );
00491 }
00492 }
00493
00494 fprintf(ioQQQ," all strong heating lines \n line wl ener/total\n");
00495 chInfo = 'h';
00496 fprintf( ioQQQ, " Sum of heat (abs) = %.2e (rel)= %.2e\n",
00497 thermal.power, thermal.power/MAX2(1e-30,total) );
00498 for( i=0; i < LineSave.nsum; i++ )
00499 {
00500 if( LineSv[i].chSumTyp == chInfo && LineSv[i].SumLine[0]/total > RATIO )
00501 {
00502 fprintf( ioQQQ, " %4.4s ", LineSv[i].chALab );
00503 prt_wl( ioQQQ, LineSv[i].wavelength );
00504 fprintf( ioQQQ, " %7.3f %c\n", LineSv[i].SumLine[0]/total, chInfo );
00505 }
00506 }
00507
00508 return;
00509 }