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

Generated for cloudy by doxygen 1.8.3.1