cloudy trunk
Loading...
Searching...
No Matches
hcmap.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/*map_do produce map of heating-cooling space for specified zone, called as result of
4 * map command */
5#include "cddefines.h"
6#include "thermal.h"
7#include "cooling.h"
8#include "called.h"
9#include "dense.h"
10#include "mole.h"
11#include "phycon.h"
12#include "trace.h"
13#include "pressure.h"
14#include "conv.h"
15#include "hcmap.h"
16#ifdef EPS
17# undef EPS
18#endif
19#define EPS 0.005
20
22
23void map_do(
24 FILE *io,
25 const char *chType)
26{
27 char chLabel[NCOLNT_LAB_LEN+1];
28 char units;
29 long int i,
30 ksav,
31 j,
32 jsav,
33 k,
34 nelem;
35 realnum wl;
36 double cfrac,
37 ch,
38 fac,
39 factor,
40 hfrac,
41 oldch,
42 ratio,
43 strhet,
44 strong,
45 t1,
46 tlowst,
47 tmax,
48 tmin,
49 torg;
50
51 DEBUG_ENTRY( "map_do()" );
52
53 t1 = phycon.te;
54 torg = phycon.te;
55 hcmap.lgMapOK = true;
56 /* flag indicating that we have computed a map */
57 hcmap.lgMapDone = true;
58
59 /* make sure pressure has been evaluated */
60 /* this sets values of pressure.PresTotlCurr */
62
63 /* print out all coolants if all else fails */
64 if( called.lgTalk )
65 {
66 fprintf( io, " Cloudy punts, Te=%10.3e HTOT=%10.3e CTOT=%10.3e nzone=%4ld\n",
67 phycon.te, thermal.htot, thermal.ctot, nzone );
68 fprintf( io, " COOLNG array is\n" );
69
70 if( thermal.ctot > 0. )
71 {
72 coolpr(io, "ZERO",1,0.,"ZERO");
73 for( i=0; i < thermal.ncltot; i++ )
74 {
75 ratio = thermal.cooling[i]/thermal.ctot;
76 if( ratio>EPS )
77 {
78 coolpr(io, (char*)thermal.chClntLab[i],thermal.collam[i],
79 ratio,"DOIT");
80 }
81 }
82
83 coolpr(io, "DONE",1,0.,"DONE");
84 fprintf( io, " Line heating array follows\n" );
85 coolpr(io, "ZERO",1,0.,"ZERO");
86
87 for( i=0; i < thermal.ncltot; i++ )
88 {
89 ratio = thermal.heatnt[i]/thermal.ctot;
90 if( ratio>EPS )
91 {
92 coolpr(io, (char*)thermal.chClntLab[i],thermal.collam[i],
93 ratio,"DOIT");
94 }
95 }
96
97 coolpr(io,"DONE",1,0.,"DONE");
98 }
99 }
100
101 /* map out te-ionization-cooling space before punching out. */
102 if( called.lgTalk )
103 {
104 fprintf( io, " map of heating, cooling, vs temp, follows.\n");
105 fprintf( io,
106 " Te Heat--------------------> Cool---------------------> dH/dT dC/DT Ne NH HII Helium \n" );
107 }
108
109 if( strcmp(chType,"punt") == 0 )
110 {
111 /* this is the original use of punt, we are punting
112 * only map within factor of two of final temperature
113 * fac will the range to either side of punted temperature */
114 fac = 1.5;
115 tmin = torg/fac;
116 tmax = torg*fac;
117
118 /* we want about 20 steps between high and low temperature
119 * default of nMapStep is 20, set with set nmaps command */
120 factor = pow(10.,log10(tmax/tmin)/(double)(hcmap.nMapStep));
121 TempChange(tmin , false);
122 }
123
124 else if( strcmp(chType," map") == 0 )
125 {
126 /* create some sort of map of heating-cooling */
127 tlowst = MAX2(hcmap.RangeMap[0],phycon.TEMP_LIMIT_LOW);
128 tmin = tlowst*0.998;
129 tmax = MIN2(hcmap.RangeMap[1],phycon.TEMP_LIMIT_HIGH)*1.002;
130
131 /* we want about nMapStep (=20) steps between high and low temperature */
132 factor = pow(10.,log10(tmax/tmin)/(double)(hcmap.nMapStep));
133 double TeNew;
134 if( thermal.lgTeHigh )
135 {
136 /* high te */
137 factor = 1./factor;
138 /* TeHighest is highest possible temperature, 1E10 */
139 TeNew = (MIN2(hcmap.RangeMap[1],phycon.TEMP_LIMIT_HIGH)/factor);
140 }
141
142 else
143 {
144 /* low te */
145 TeNew = (tlowst/factor);
146 }
147 TempChange(TeNew , false);
148 }
149
150 else
151 {
152 /* don't know what to do */
153 fprintf( ioQQQ, " PUNT called with insane argument,=%4.4s\n",
154 chType );
156 }
157
158 /* now allocate space for te, c, h vectors in map, if not already done */
159 if( hcmap.nMapAlloc==0 )
160 {
161 /* space not allocated, do so now */
162 hcmap.nMapAlloc = hcmap.nMapStep+4;
163
164 /* now make the space */
165 hcmap.temap = (realnum *)MALLOC( sizeof(realnum)*(size_t)(hcmap.nMapStep+4) );
166 if( hcmap.temap == NULL )
167 {
168 printf( " not enough memory to allocate hcmap.temap in map_do\n" );
170 }
171 hcmap.cmap = (realnum *)MALLOC( sizeof(realnum)*(size_t)(hcmap.nMapStep+4) );
172 if( hcmap.cmap == NULL )
173 {
174 printf( " not enough memory to allocate hcmap.cmap in map_do\n" );
176 }
177 hcmap.hmap = (realnum *)MALLOC( sizeof(realnum)*(size_t)(hcmap.nMapStep+4) );
178 if( hcmap.hmap == NULL )
179 {
180 printf( " not enough memory to allocate hcmap.hmap in map_do\n" );
182 }
183 }
184
185 thermal.lgCNegChk = false;
186 hcmap.nmap = 0;
187 oldch = 0.;
188 TempChange(phycon.te *factor , true);
189 if( trace.nTrConvg )
190 fprintf(ioQQQ, " MAP called temp range %.4e %.4e in %li stops ===============================================\n",
191 tmin,
192 tmax,
193 hcmap.nmap);
194
195 while( (double)phycon.te < tmax*0.999 && (double)phycon.te > tmin*1.001 )
196 {
197 /* this sets values of pressure.PresTotlCurr */
199
200 /* must reset upper and lower bounds for ionization distributions */
201 /* fix number of stages of ionization */
202 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
203 {
204 if( dense.lgElmtOn[nelem] )
205 {
206 dense.IonLow[nelem] = 0;
207 /*
208 * IonHigh[n] is the highest stage of ionization present
209 * the IonHigh array index is on the C scake, so [0] is hydrogen
210 * the value is also on the C scale, so element [nelem] can range
211 * from 0 to nelem+1
212 */
213 dense.IonHigh[nelem] = nelem + 1;
214 }
215 else
216 {
217 /* this element is turned off, make stages impossible */
218 dense.IonLow[nelem] = -1;
219 dense.IonHigh[nelem] = -1;
220 }
221 }
222
223 /* this turns on constant reevaluation of everything */
224 conv.lgSearch = true;
225
226 if( trace.nTrConvg )
227 fprintf(ioQQQ, " MAP new temp %.4e ===============================================\n",
228 phycon.te );
229
230 /* this counts how many times ionize is called in this model after startr,
231 * and is flag used by ionize to understand it is being called the first time*/
232 conv.nTotalIoniz = 0;
233
234 /* now get ionization solution for this temperature */
235 if( ConvEdenIoniz() )
236 lgAbort = true;
237
238 /* save results for later prints */
239 hcmap.temap[hcmap.nmap] = (realnum)phycon.te;
240 hcmap.cmap[hcmap.nmap] = (realnum)thermal.ctot;
241 hcmap.hmap[hcmap.nmap] = (realnum)thermal.htot;
242
243 wl = 0.f;
244 strong = 0.;
245
246 for( j=0; j < thermal.ncltot; j++ )
247 {
248 if( thermal.cooling[j] > strong )
249 {
250 strcpy( chLabel, thermal.chClntLab[j] );
251 strong = thermal.cooling[j];
252 wl = thermal.collam[j];
253 }
254 }
255
256 cfrac = strong/thermal.ctot;
257 strhet = 0.;
258 /* these will be reset in following loop*/
259 ksav = -INT_MAX;
260 jsav = -INT_MAX;
261
262 for( k=0; k < LIMELM; k++ )
263 {
264 for( j=0; j < LIMELM; j++ )
265 {
266 if( thermal.heating[k][j] > strhet )
267 {
268 strhet = thermal.heating[k][j];
269 jsav = j;
270 ksav = k;
271 }
272 }
273 }
274
275 ch = thermal.ctot - thermal.htot;
276 /* use ratio to check for change of sign since product
277 * can underflow at low densities */
278 if( oldch/ch < 0. && called.lgTalk )
279 {
280 fprintf( io, " ----------------------------------------------- Probable thermal solution here. --------------------------------------------\n" );
281 }
282
283 oldch = ch;
284 hfrac = strhet/thermal.htot;
285 if( called.lgTalk )
286 {
287 /* convert to micros if IR transition */
288 if( wl < 100000.f )
289 {
290 units = 'A';
291 }
292
293 else
294 {
295 wl /= 10000.f;
296 units = 'm';
297 }
298
299 if( trace.lgTrace )
300 {
301 fprintf( io, "TRACE: te, htot, ctot%11.3e%11.3e%11.3e\n",
302 phycon.te, thermal.htot, thermal.ctot );
303 }
304
305 /*fprintf( io,
306 "%10.4e%11.4e%4ld%4ld%6.3f%11.4e %4.4s %4ld%c%6.3f%10.2e%11.4e%11.4e%6.2f%6.2f%6.2f%6.2f\n",;*/
307 fprintf(io, PrintEfmt("%11.4e", phycon.te ) );
308 fprintf(io, PrintEfmt("%11.4e", thermal.htot ) );
309 fprintf(io," [%2ld][%2ld]%6.3f",
310 ksav, jsav,
311 hfrac);
312 fprintf(io, PrintEfmt("%11.4e", thermal.ctot ) );
313 fprintf(io," %s %.1f%c%6.3f",
314 chLabel ,
315 wl,
316 units,
317 cfrac );
318 fprintf(io, PrintEfmt(" %10.2e", thermal.dHeatdT ) );
319 fprintf(io, PrintEfmt("%11.2e", thermal.dCooldT ) );
320 fprintf(io, PrintEfmt("%11.4e", dense.eden ) );
321 fprintf(io, PrintEfmt("%11.4e", dense.gas_phase[ipHYDROGEN] ) );
322 if( dense.lgElmtOn[ipHELIUM] )
323 {
324 fprintf(io,"%6.2f%6.2f%6.2f%6.2f\n",
325 log10(MAX2(1e-9,dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN])),
326 log10(MAX2(1e-9,dense.xIonDense[ipHELIUM][0]/dense.gas_phase[ipHELIUM])),
327 log10(MAX2(1e-9,dense.xIonDense[ipHELIUM][1]/dense.gas_phase[ipHELIUM])),
328 log10(MAX2(1e-9,dense.xIonDense[ipHELIUM][2]/dense.gas_phase[ipHELIUM])) );
329 }
330 fflush(io);
331 }
332
333 TempChange(phycon.te*factor , true);
334 /* increment nmap but do not exceed nMapAlloc */
335 hcmap.nmap = MIN2(hcmap.nMapAlloc,hcmap.nmap+1);
336
337 {
338 enum {DEBUG_LOC=false};
339 if( DEBUG_LOC )
340 {
341 static int kount = 0;
342 factor = 1.;
343 TempChange(8674900. , true);
344 ++kount;
345 if( kount >=100 )
346 {
347 fprintf(ioQQQ," exiting in map_do\n");
348 break;
349 }
350 }
351 }
352 }
353
354 /* now check whether sharp inflections occurred, and also find the biggest jump
355 * in the heating and cooling */
356 hcmap.lgMapOK = true;
357 /* >>chng 02 mar 04, lower bound had been 1, so [i-2] below was negative */
358 for( i=2; i< hcmap.nmap-2; ++i )
359 {
360 realnum s1,s2,s3;/* the three slopes we will use */
361 s1 = hcmap.cmap[i-2] - hcmap.cmap[i-1];
362 s2 = hcmap.cmap[i-1] - hcmap.cmap[i];
363 s3 = hcmap.cmap[i] - hcmap.cmap[i+1];
364 if( s1*s3 > 0. && s2*s3 < 0. )
365 {
366 /* of the three points, the outer had the same slope
367 * (their product was positive) but there was an inflection
368 * between them (the negative product). The data chain looked like
369 * 2 4
370 * 1 3 or vice versa, either case is wrong,
371 * with the logic in the above test, the problem point will aways be s2 */
372 fprintf( io,
373 " cooling curve had double inflection at T=%.2e. ",
374 hcmap.temap[i]);
375 fprintf( io, " Slopes were %.2e %.2e %.2e", s1, s2, s3);
376 if( fabs(s2)/hcmap.cmap[i] > 0.05 )
377 {
378 fprintf( io,
379 " error large, (rel slope of %.2e).\n",
380 s2 / hcmap.cmap[i]);
381 hcmap.lgMapOK = false;
382 }
383 else
384 {
385 fprintf( io,
386 " error is small, (rel slope of %.2e).\n",
387 s2 / hcmap.cmap[i]);
388 }
389 }
390
391 s1 = hcmap.hmap[i-2] - hcmap.hmap[i-1];
392 s2 = hcmap.hmap[i-1] - hcmap.hmap[i];
393 s3 = hcmap.hmap[i] - hcmap.hmap[i+1];
394 if( s1*s3 > 0. && s2*s3 < 0. )
395 {
396 /* of the three points, the outer had the same slope
397 * (their product was positive) but there was an inflection
398 * between them (the negative product). The data chain looked like
399 * 2 4
400 * 1 3 or vice versa, either case is wrong */
401 fprintf( io,
402 " heating curve had double inflection at T=%.2e.\n",
403 hcmap.temap[i] );
404 hcmap.lgMapOK = false;
405 }
406 }
407
408 thermal.lgCNegChk = true;
409 TempChange(t1 , false);
410 return;
411}
t_called called
Definition called.cpp:5
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
bool lgAbort
Definition cddefines.cpp:10
#define PrintEfmt(F, V)
Definition cddefines.h:1472
#define MIN2
Definition cddefines.h:761
#define MALLOC(exp)
Definition cddefines.h:501
const int LIMELM
Definition cddefines.h:258
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int ipHELIUM
Definition cddefines.h:306
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
t_conv conv
Definition conv.cpp:5
int ConvEdenIoniz(void)
#define EPS
void coolpr(FILE *io, const char *chLabel, realnum lambda, double ratio, const char *chJOB)
Definition cool_pr.cpp:9
t_dense dense
Definition dense.cpp:24
t_hcmap hcmap
Definition hcmap.cpp:21
void map_do(FILE *io, const char *chType)
Definition hcmap.cpp:23
t_phycon phycon
Definition phycon.cpp:6
void PresTotCurrent(void)
void TempChange(double TempNew, bool lgForceUpdate)
t_thermal thermal
Definition thermal.cpp:5
#define NCOLNT_LAB_LEN
Definition thermal.h:91
t_trace trace
Definition trace.cpp:5