cloudy trunk
Loading...
Searching...
No Matches
energy.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3
4#include "cddefines.h"
5#include "energy.h"
6#include "continuum.h"
7#include "rfield.h"
8#include "doppvel.h"
9#include "geometry.h"
10#include "hextra.h"
11#include "ipoint.h"
12#include "lines.h"
13#include "lines_service.h"
14#include "opacity.h"
15#include "physconst.h"
16#include "prt.h"
17#include "radius.h"
18#include "secondaries.h"
19#include "struc.h"
20#include "thermal.h"
21#include "warnings.h"
22#include "wind.h"
23#include "dynamics.h"
24
25static const char *ENERGY_RYD = "Ryd";
26static const char *ENERGY_ERG = "erg";
27static const char *ENERGY_EV = "eV";
28static const char *ENERGY_KEV = "keV";
29static const char *ENERGY_MEV = "MeV";
30static const char *ENERGY_WN = "cm^-1";
31static const char *ENERGY_CM = "cm";
32static const char *ENERGY_MM = "mm";
33static const char *ENERGY_MICRON = "um";
34static const char *ENERGY_NM = "nm";
35static const char *ENERGY_A = "A";
36static const char *ENERGY_HZ = "Hz";
37static const char *ENERGY_KHZ = "kHz";
38static const char *ENERGY_MHZ = "MHz";
39static const char *ENERGY_GHZ = "GHz";
40static const char *ENERGY_K = "K";
41
42inline bool isSameUnit(const char *unit, const char *ref)
43{
44 return strcmp(unit,ref) == 0;
45}
46
47const char *StandardEnergyUnit(const char *chCard)
48{
49 DEBUG_ENTRY( "StandardEnergyUnit()" );
50
51 // use " MIC" so that it cannot pick up on "W/SQM/MICRON" in flux units
52 if( nMatch(" MIC",chCard) )
53 {
54 /* micron */
55 return ENERGY_MICRON;
56 }
57 else if( nMatch(" EV ",chCard) )
58 {
59 /* eV */
60 return ENERGY_EV;
61 }
62 else if( nMatch(" KEV",chCard) )
63 {
64 /* keV */
65 return ENERGY_KEV;
66 }
67 else if( nMatch(" MEV",chCard) )
68 {
69 /* MeV */
70 return ENERGY_MEV;
71 }
72 else if( nMatch("WAVE",chCard) )
73 {
74 /* wavenumbers */
75 return ENERGY_WN;
76 }
77 else if( nMatch("CENT",chCard) || nMatch(" CM ",chCard) )
78 {
79 /* centimeter */
80 return ENERGY_CM;
81 }
82 else if( nMatch(" MM ",chCard) )
83 {
84 /* millimeter */
85 return ENERGY_MM;
86 }
87 else if( nMatch(" NM ",chCard) )
88 {
89 /* nanometer */
90 return ENERGY_NM;
91 }
92 else if( nMatch("ANGS",chCard) )
93 {
94 /* angstrom */
95 return ENERGY_A;
96 }
97 else if( nMatch(" HZ ",chCard) )
98 {
99 /* hertz */
100 return ENERGY_HZ;
101 }
102 else if( nMatch(" KHZ",chCard) )
103 {
104 /* kilohertz */
105 return ENERGY_KHZ;
106 }
107 else if( nMatch(" MHZ",chCard) )
108 {
109 /* megahertz */
110 return ENERGY_MHZ;
111 }
112 else if( nMatch(" GHZ",chCard) )
113 {
114 /* gigahertz */
115 return ENERGY_GHZ;
116 }
117 else if( nMatch("KELV",chCard) || nMatch(" K ",chCard) )
118 {
119 /* kelvin */
120 return ENERGY_K;
121 }
122 else if( nMatch(" RYD",chCard) )
123 {
124 /* RydbERG must come before ERG to avoid false trip */
125 return ENERGY_RYD;
126 }
127 // use "ERG " so that it cannot pick up on "ERG/S" in flux units
128 else if( nMatch(" ERG ",chCard) )
129 {
130 /* erg */
131 return ENERGY_ERG;
132 }
133 else
134 {
135 fprintf( ioQQQ, " No energy / wavelength unit was recognized on this line:\n %s\n\n", chCard );
136 fprintf( ioQQQ, " See Hazy for details.\n" );
138 }
139}
140
141double Energy::get(const char *unit) const
142{
143 DEBUG_ENTRY( "Energy::get()" );
144
145 /* internal unit is Ryd */
146 if( isSameUnit(unit,ENERGY_RYD) )
147 {
148 return Ryd();
149 }
150 else if( isSameUnit(unit,ENERGY_ERG) )
151 {
152 return Erg();
153 }
154 else if( isSameUnit(unit,ENERGY_EV) )
155 {
156 return eV();
157 }
158 else if( isSameUnit(unit,ENERGY_KEV) )
159 {
160 return keV();
161 }
162 else if( isSameUnit(unit,ENERGY_MEV) )
163 {
164 return MeV();
165 }
166 else if( isSameUnit(unit,ENERGY_WN) )
167 {
168 return WN();
169 }
170 else if( isSameUnit(unit,ENERGY_CM) )
171 {
172 return cm();
173 }
174 else if( isSameUnit(unit,ENERGY_MM) )
175 {
176 return mm();
177 }
178 else if( isSameUnit(unit,ENERGY_MICRON) )
179 {
180 return micron();
181 }
182 else if( isSameUnit(unit,ENERGY_NM) )
183 {
184 return nm();
185 }
186 else if( isSameUnit(unit,ENERGY_A) )
187 {
188 return Angstrom();
189 }
190 else if( isSameUnit(unit,ENERGY_HZ) )
191 {
192 return Hz();
193 }
194 else if( isSameUnit(unit,ENERGY_KHZ) )
195 {
196 return kHz();
197 }
198 else if( isSameUnit(unit,ENERGY_MHZ) )
199 {
200 return MHz();
201 }
202 else if( isSameUnit(unit,ENERGY_GHZ) )
203 {
204 return GHz();
205 }
206 else if( isSameUnit(unit,ENERGY_K) )
207 {
208 return K();
209 }
210 else
211 {
212 fprintf( ioQQQ, " insane units in Energy::get: \"%s\"\n", unit );
214 }
215}
216
217void Energy::set(double value, const char *unit)
218{
219 DEBUG_ENTRY( "Energy::set()" );
220
221 /* internal unit is Ryd */
222 if( isSameUnit(unit,ENERGY_RYD) )
223 {
224 set(value);
225 }
226 else if( isSameUnit(unit,ENERGY_ERG) )
227 {
228 set(value/EN1RYD);
229 }
230 else if( isSameUnit(unit,ENERGY_MEV) )
231 {
232 set(1e6*value/EVRYD);
233 }
234 else if( isSameUnit(unit,ENERGY_KEV) )
235 {
236 set(1e3*value/EVRYD);
237 }
238 else if( isSameUnit(unit,ENERGY_EV) )
239 {
240 set(value/EVRYD);
241 }
242 else if( isSameUnit(unit,ENERGY_WN) )
243 {
244 set(value/RYD_INF);
245 }
246 else if( isSameUnit(unit,ENERGY_A) )
247 {
248 set(RYDLAM/value);
249 }
250 else if( isSameUnit(unit,ENERGY_NM) )
251 {
252 set(RYDLAM/(1.e1*value));
253 }
254 else if( isSameUnit(unit,ENERGY_MICRON) )
255 {
256 set(RYDLAM/(1.e4*value));
257 }
258 else if( isSameUnit(unit,ENERGY_MM) )
259 {
260 set(RYDLAM/(1.e7*value));
261 }
262 else if( isSameUnit(unit,ENERGY_CM) )
263 {
264 set(RYDLAM/(1.e8*value));
265 }
266 else if( isSameUnit(unit,ENERGY_HZ) )
267 {
268 set(value/FR1RYD);
269 }
270 else if( isSameUnit(unit,ENERGY_KHZ) )
271 {
272 set(1e3*value/FR1RYD);
273 }
274 else if( isSameUnit(unit,ENERGY_MHZ) )
275 {
276 set(1e6*value/FR1RYD);
277 }
278 else if( isSameUnit(unit,ENERGY_GHZ) )
279 {
280 set(1e9*value/FR1RYD);
281 }
282 else if( isSameUnit(unit,ENERGY_K) )
283 {
284 set(value/TE1RYD);
285 }
286 else
287 {
288 fprintf( ioQQQ, " insane units in Energy::set: \"%s\"\n", unit );
290 }
291}
292
294{
295 DEBUG_ENTRY( "EnergyEntry::p_set_ip()" );
296
297 double energy = Ryd();
298 if( energy < rfield.emm || energy > rfield.egamry )
299 {
300 fprintf( ioQQQ, " The energy %g Ryd is not within the default Cloudy range\n", energy );
302 }
303 p_ip = ipoint(energy) - 1;
304 ASSERT( p_ip >= 0 );
305}
306
307STATIC double sum_radiation( void );
308
309/*badprt print out coolants if energy not conserved */
310STATIC void badprt(
311 /* total is total luminosity available in radiation */
312 double total);
313
315{
316 double outin,
317 outout,
318 outtot,
319 sum_coolants,
320 sum_recom_lines,
321 flow_energy;
322
323 char chLine[INPUT_LINE_LENGTH];
324
325 DEBUG_ENTRY( "lgConserveEnergy()" );
326
327 /* check whether energy is conserved
328 * following is outward continuum */
329 outsum(&outtot,&outin,&outout);
330 /* sum_recom_lines is the sum of all recombination line energy */
331 sum_recom_lines = totlin('r');
332 if( sum_recom_lines == 0. )
333 {
334 sprintf( chLine, " !Total recombination line energy is 0." );
335 bangin(chLine);
336 }
337
338 /* sum_coolants is the sum of all coolants */
339 sum_coolants = totlin('c');
340 if( sum_coolants == 0. )
341 {
342 sprintf( chLine, " !Total cooling is zero." );
343 bangin(chLine);
344 }
345
346 if( !(wind.lgBallistic() || wind.lgStatic()) )
347 {
348 /* approximate energy density coming out in wind
349 * should ideally include flow into grid and internal energies */
350 flow_energy = (2.5*struc.GasPressure[0]+0.5*struc.DenMass[0]*wind.windv0*wind.windv0)
351 *(-wind.windv0);
352 }
353 else
354 {
355 flow_energy = 0.;
356 }
357
358 if( 0 && geometry.lgSphere && (!thermal.lgTemperatureConstant) && (!dynamics.lgTimeDependentStatic) &&
359 (hextra.cryden == 0.) && ((hextra.TurbHeat+DoppVel.DispScale) == 0.) && !secondaries.lgCSetOn )
360 {
361 double sum_rad = sum_radiation();
362
363 if( fabs( 1. - sum_rad/(continuum.TotalLumin*POW2(radius.rinner)) ) > 0.01 )
364 fprintf( ioQQQ, "PROBLEM DISASTER lgConserveEnergy failed %li\t%li\t%e\t%e\t%e\n",
365 iteration, nzone, sum_rad, continuum.TotalLumin*POW2(radius.rinner),
366 sum_rad/(continuum.TotalLumin*POW2(radius.rinner)) );
367 }
368
369 /* TotalLumin is total incident photon luminosity, evaluated in setup
370 * sum_coolants is evaluated here, is sum of all coolants */
371 /* this test is meaningless when the APERTURE command is in effect since
372 * sum_coolants and sum_recom_lines react to the APERTURE command, while
373 * continuum.TotalLumin does not */
374 if( !dynamics.lgTimeDependentStatic && (sum_coolants + sum_recom_lines + flow_energy) > continuum.TotalLumin*geometry.covgeo &&
375 !thermal.lgTemperatureConstant && geometry.iEmissPower == 2 )
376 {
377 if( (hextra.cryden == 0.) && ((hextra.TurbHeat+DoppVel.DispScale) == 0.) && !secondaries.lgCSetOn )
378 {
379
380 sprintf( chLine,
381 " W-Radiated luminosity (cool+rec+wind=%.2e+%.2e+%.2e) is greater than that in incident radiation (TotalLumin=%8.2e). Power radiated was %8.2e",
382 sum_coolants, sum_recom_lines, flow_energy , continuum.TotalLumin, thermal.power );
383 warnin(chLine);
384 /* write same thing directly to output (above will be sorted later) */
385 fprintf( ioQQQ, "\n\n DISASTER This calculation DID NOT CONSERVE ENERGY!\n\n\n" );
386
387 if( !continuum.lgCheckEnergyEveryZone )
388 fprintf( ioQQQ, "Rerun with *set check energy every zone* command to do energy check for each zone.\n\n");
389
390 lgAbort = true;
391
392 /* the case b command can cause this problem - say so if case b was set */
393 if( opac.lgCaseB )
394 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" );
395
396 /* print out significant heating and cooling */
397 badprt(continuum.TotalLumin);
398
399 sprintf( chLine, " W-Something is really wrong: the ratio of radiated to incident luminosity is %.2e.",
400 (sum_coolants + sum_recom_lines)/continuum.TotalLumin );
401 warnin(chLine);
402
403 /* this can be caused by the force te command */
404 if( thermal.ConstTemp > 0. )
405 {
406 fprintf( ioQQQ, " This may have been caused by the FORCE TE command.\n" );
407 fprintf( ioQQQ, " Remove it and run again.\n" );
408 }
409 else
410 return false;
411 }
412 }
413
414 return true;
415}
416
417STATIC double sum_radiation( void )
418{
419 double flxref = 0., flxatt = 0., conem = 0.;
420 long nEmType = 0;
421 double sum = 0.;
422
423 const realnum *trans_coef_total=rfield.getCoarseTransCoef();
424 for( long j=0; j<rfield.nflux; j++ )
425 {
426 /* the reflected continuum */
427 flxref += (rfield.ConRefIncid[nEmType][j] + rfield.ConEmitReflec[nEmType][j] + rfield.reflin[nEmType][j]) * rfield.anu[j];
428 ASSERT( flxref >= 0. );
429
430 /* the attenuated incident continuum */
431 flxatt += rfield.flux[nEmType][j]*radius.r1r0sq*trans_coef_total[j] * rfield.anu[j];
432 ASSERT( flxatt >= 0. );
433
434 /* the outward emitted continuum */
435 conem += (rfield.ConEmitOut[nEmType][j] + rfield.outlin[nEmType][j] + rfield.outlin_noplot[j])*radius.r1r0sq*geometry.covgeo * rfield.anu[j];
436 ASSERT( conem >= 0. );
437 }
438
439 sum = (flxref + flxatt + conem) * EN1RYD * POW2(radius.rinner);
440 ASSERT( sum >= 0. );
441
442 return sum;
443}
444
445/*badprt print out coolants if energy not conserved */
447 /* total is total luminosity available in radiation */
448 double total)
449{
450 /* following is smallest ratio to print */
451 static const double RATIO = 0.02;
452 char chInfo;
453 long int i;
454 realnum sum_coolants,
455 sum_recom_lines;
456
457 DEBUG_ENTRY( "badprt()" );
458
459 fprintf( ioQQQ, " badprt: all entries with greater than%6.2f%% of incident continuum follow.\n",
460 RATIO*100. );
461 fprintf( ioQQQ, " badprt: Intensities are relative to total energy in incident continuum.\n" );
462
463 /* now find sum of recombination lines */
464 chInfo = 'r';
465 sum_recom_lines = (realnum)totlin('r');
466 fprintf( ioQQQ,
467 " Sum of energy in recombination lines is %.2e relative to total incident radiation is %.2e\n",
468 sum_recom_lines,
469 sum_recom_lines/MAX2(1e-30,total) );
470
471 fprintf(ioQQQ," all strong information lines \n line wl ener/total\n");
472 /* now print all strong lines */
473 for( i=0; i < LineSave.nsum; i++ )
474 {
475 if( LineSv[i].chSumTyp == chInfo && LineSv[i].SumLine[0]/total > RATIO )
476 {
477 fprintf( ioQQQ, " %4.4s ", LineSv[i].chALab );
479 fprintf( ioQQQ, " %7.3f %c\n", LineSv[i].SumLine[0]/total, chInfo );
480 }
481 }
482
483 fprintf(ioQQQ," all strong cooling lines \n line wl ener/total\n");
484 chInfo = 'c';
485 sum_coolants = (realnum)totlin('c');
486 fprintf( ioQQQ, " Sum of coolants (abs) = %.2e (rel)= %.2e\n",
487 sum_coolants, sum_coolants/MAX2(1e-30,total) );
488 for( i=0; i < LineSave.nsum; i++ )
489 {
490 if( LineSv[i].chSumTyp == chInfo && LineSv[i].SumLine[0]/total > RATIO )
491 {
492 fprintf( ioQQQ, " %4.4s ", LineSv[i].chALab );
494 fprintf( ioQQQ, " %7.3f %c\n", LineSv[i].SumLine[0]/total, chInfo );
495 }
496 }
497
498 fprintf(ioQQQ," all strong heating lines \n line wl ener/total\n");
499 chInfo = 'h';
500 fprintf( ioQQQ, " Sum of heat (abs) = %.2e (rel)= %.2e\n",
501 thermal.power, thermal.power/MAX2(1e-30,total) );
502 for( i=0; i < LineSave.nsum; i++ )
503 {
504 if( LineSv[i].chSumTyp == chInfo && LineSv[i].SumLine[0]/total > RATIO )
505 {
506 fprintf( ioQQQ, " %4.4s ", LineSv[i].chALab );
508 fprintf( ioQQQ, " %7.3f %c\n", LineSv[i].SumLine[0]/total, chInfo );
509 }
510 }
511
512 return;
513}
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
bool lgAbort
Definition cddefines.cpp:10
long int iteration
Definition cddefines.cpp:16
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
long nMatch(const char *chKey, const char *chCard)
Definition service.cpp:451
#define POW2
Definition cddefines.h:929
const int INPUT_LINE_LENGTH
Definition cddefines.h:254
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
LinSv * LineSv
Definition cdinit.cpp:70
void p_set_ip()
Definition energy.cpp:293
long p_ip
Definition energy.h:100
double GHz() const
Definition energy.h:66
double mm() const
Definition energy.h:82
double MHz() const
Definition energy.h:62
double MeV() const
Definition energy.h:50
double get(const char *unit) const
Definition energy.cpp:141
double Ryd() const
Definition energy.h:26
double Angstrom() const
Definition energy.h:70
double nm() const
Definition energy.h:74
double WN() const
Definition energy.h:38
double cm() const
Definition energy.h:86
double Hz() const
Definition energy.h:54
double kHz() const
Definition energy.h:58
double keV() const
Definition energy.h:46
double micron() const
Definition energy.h:78
double K() const
Definition energy.h:30
double Erg() const
Definition energy.h:34
void set(double energy)
Definition energy.h:19
double eV() const
Definition energy.h:42
void outsum(double *outtot, double *outin, double *outout)
long ipoint(double energy_ryd)
t_continuum continuum
Definition continuum.cpp:5
t_DoppVel DoppVel
Definition doppvel.cpp:5
t_dynamics dynamics
Definition dynamics.cpp:44
static const char * ENERGY_CM
Definition energy.cpp:31
static const char * ENERGY_WN
Definition energy.cpp:30
static const char * ENERGY_KHZ
Definition energy.cpp:37
bool lgConserveEnergy()
Definition energy.cpp:314
static const char * ENERGY_RYD
Definition energy.cpp:25
STATIC double sum_radiation(void)
Definition energy.cpp:417
static const char * ENERGY_MEV
Definition energy.cpp:29
static const char * ENERGY_A
Definition energy.cpp:35
static const char * ENERGY_K
Definition energy.cpp:40
static const char * ENERGY_GHZ
Definition energy.cpp:39
STATIC void badprt(double total)
Definition energy.cpp:446
static const char * ENERGY_MICRON
Definition energy.cpp:33
const char * StandardEnergyUnit(const char *chCard)
Definition energy.cpp:47
static const char * ENERGY_KEV
Definition energy.cpp:28
static const char * ENERGY_ERG
Definition energy.cpp:26
static const char * ENERGY_MHZ
Definition energy.cpp:38
static const char * ENERGY_EV
Definition energy.cpp:27
static const char * ENERGY_NM
Definition energy.cpp:34
static const char * ENERGY_HZ
Definition energy.cpp:36
bool isSameUnit(const char *unit, const char *ref)
Definition energy.cpp:42
static const char * ENERGY_MM
Definition energy.cpp:32
t_geometry geometry
Definition geometry.cpp:5
t_hextra hextra
Definition hextra.cpp:5
t_LineSave LineSave
Definition lines.cpp:5
double totlin(int chInfo)
static realnum * wavelength
t_opac opac
Definition opacity.cpp:5
UNUSED const double FR1RYD
Definition physconst.h:195
UNUSED const double RYD_INF
Definition physconst.h:115
UNUSED const double EN1RYD
Definition physconst.h:179
UNUSED const double EVRYD
Definition physconst.h:189
UNUSED const double TE1RYD
Definition physconst.h:183
UNUSED const double RYDLAM
Definition physconst.h:176
void prt_wl(FILE *ioOUT, realnum wl)
Definition prt.cpp:13
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
t_secondaries secondaries
t_struc struc
Definition struc.cpp:6
t_thermal thermal
Definition thermal.cpp:5
void warnin(char *chLine)
Definition warnings.cpp:27
void bangin(char *chLine)
Definition warnings.cpp:73
Wind wind
Definition wind.cpp:5