cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
plot.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 /*plot master routine to generate some sort of plot */
4 #include "cddefines.h"
5 #include "iterations.h"
6 #include "called.h"
7 #define IHI 59
8 #define IWID 121
9 #include "input.h"
10 #include "rfield.h"
11 #include "trace.h"
12 #include "radius.h"
13 #include "geometry.h"
14 #include "opacity.h"
15 #include "dense.h"
16 #include "hcmap.h"
17 #include "plot.h"
18 
19 /*pltcon generate plot of continuum array */
20 STATIC void pltcon(long int np,
21  const char *chCall);
22 
23 /*pltmap generate plot of heating and cooling map */
24 STATIC void pltmap(long int np,
25  const char *chCall);
26 
27 /*pltopc generate plot of local gas opacity */
28 STATIC void pltopc(long int np,
29  const char *chCall);
30 
31 /* this is the base routine that actually makes the plots, called by above */
32 STATIC void pltr(realnum[],realnum[],long,double,double,double,double,
33  char,char*,long,bool);
34 
35 
36 void plot(const char *chCall)
37 {
38  long int np;
39 
40  DEBUG_ENTRY( "plot()" );
41 
42  /* return if this is not the last iteration, or a plot not required,
43  * or we are not speaking */
44  if( !plotCom.lgPlotON || !called.lgTalk )
45  {
46  return;
47  }
48 
49  if( !iterations.lgLastIt && (strcmp(chCall,"FIRST") != 0) )
50  {
51  return;
52  }
53 
54  /* loop over all the requested plots */
55  for( np=0; np < plotCom.nplot; np++ )
56  {
57  /* series of tests to determine which type of plot we will do */
58  if( strcmp(plotCom.chPType[np]," MAP") == 0 )
59  {
60  /* thermal map */
61  pltmap(np,chCall);
62  }
63  else if( strcmp(plotCom.chPType[np] ,"CONT") == 0 ||
64  strcmp(plotCom.chPType[np] ,"CRAW") == 0 ||
65  strcmp(plotCom.chPType[np] ,"DIFF") == 0 ||
66  strcmp(plotCom.chPType[np] ,"REFL") == 0 ||
67  strcmp(plotCom.chPType[np] ,"EMIT") == 0 ||
68  strcmp(plotCom.chPType[np] ,"CPHT") == 0 ||
69  strcmp(plotCom.chPType[np] ,"OUTW") == 0 )
70  {
71  /* this is a contiuum plot of some kind */
72  pltcon(np,chCall);
73  }
74 
75  else if(
76  strcmp(plotCom.chPType[np] ,"OPAA") == 0 ||
77  strcmp(plotCom.chPType[np] ,"OPAS") == 0 ||
78  strcmp(plotCom.chPType[np] ,"OPAT") == 0 )
79  {
80  /* absorption, scattering, or total opacity */
81  pltopc(np,chCall);
82  }
83  else
84  {
85  fprintf( ioQQQ, " PLOT type=%4.4s not known. STOP\n",
86  plotCom.chPType[np] );
87  cdEXIT(EXIT_FAILURE);
88  }
89  }
90 
91  return;
92 }
93 
94 /*pltcon generate plot of continuum array */
95 
97  long int np,
98  const char *chCall)
99 {
100  char chSymPlt2[3],
101  chXtitle[23];
102  char chSym,
103  chSymPlt1;
104  long int i;
105  double contin,
106  ymin2;
107  static double xmax,
108  xmin,
109  ymax,
110  ymin;
111  static realnum *y/*[rfield.nupper]*/,
112  *y2/*[rfield.nupper]*/;
113 
114  DEBUG_ENTRY( "pltcon()" );
115 
116  if( strcmp(chCall,"FIRST") == 0 )
117  {
118  return;
119  }
120 
121  xmin = rfield.anulog[0];
122  xmin = MAX2((double)plotCom.pltxmn[np],xmin);
123  xmax = rfield.anulog[rfield.nflux-1];
124  xmax = MIN2(xmax,(double)plotCom.pltxmx[np]);
125 
126  if( plotCom.lgPltTrace[np] )
127  {
128  fprintf( ioQQQ, " XMIN, XMAX=%12.4e%12.4e NFLUX=%4ld\n",
129  xmin, xmax, rfield.nflux );
130  }
131 
132  if( strcmp(plotCom.chPType[np],"REFL") == 0 && geometry.lgSphere )
133  {
134  fprintf( ioQQQ, " Reflected continuum not computed when SPHERE set.\n" );
135  return;
136  }
137 
138  y = (realnum*)MALLOC((size_t)rfield.nupper*sizeof(realnum) );
139  y2 = (realnum*)MALLOC((size_t)rfield.nupper*sizeof(realnum) );
140 
141  /* these will be the default symbols for first and second plot */
142  chSymPlt1 = '.';
143  strcpy( chSymPlt2, "o " );
144  ymin = FLT_MAX;
145  ymin2 = FLT_MAX;
146  ymax = -FLT_MAX;
147  for( i=0; i < rfield.nflux; i++ )
148  {
149  if( (double)rfield.anulog[i] > xmin && (double)rfield.anulog[i] < xmax )
150  {
151  if( strcmp(plotCom.chPType[np],"CONT") == 0 )
152  {
153  y[i] = (realnum)log10(MAX2(rfield.flux_total_incident[i]/rfield.widflx[i]*
154  rfield.anu2[i],1e-37));
155  /* >>chng 01 jul 13, add rfield.ConEmitReflec[i] */
157  y2[i] = (realnum)MAX2((contin/rfield.widflx[i]+(rfield.outlin[i]+rfield.outlin_noplot[i])/
159  rfield.anu2[i]*radius.r1r0sq,1e-37);
160  y2[i] = (realnum)log10(y2[i]);
161  }
162  else if( strcmp(plotCom.chPType[np],"CPHT") == 0 )
163  {
164  /* plot continuum as photons */
165  y[i] = (realnum)log10(MAX2(rfield.flux_total_incident[i]/rfield.widflx[i],
166  1e-37));
167  contin = rfield.flux[i] + rfield.ConEmitOut[i]*geometry.covgeo/rfield.widflx[i];
168  y2[i] = (realnum)MAX2((contin+(rfield.outlin[i]+rfield.outlin[i])/rfield.anu[i]*
169  geometry.covgeo)* radius.r1r0sq,1e-37);
170  y2[i] = (realnum)log10(y2[i]);
171  }
172  else if( strcmp(plotCom.chPType[np],"REFL") == 0 )
173  {
174  /* plot "reflected" continuum from last zone only */
175  y[i] = (realnum)log10(MAX2((rfield.ConEmitReflec[i]/rfield.widflx[i]+
176  rfield.reflin[i])*rfield.anu2[i],1e-37));
177  y2[i] = y[i];
178  }
179  else if( strcmp(plotCom.chPType[np],"EMIT") == 0 )
180  {
181  /* plot "emitted" continuum from both sides of cloud */
182  y[i] = (realnum)log10(MAX2(
184  rfield.widflx[i]+
186  rfield.anu2[i],1e-37));
187  y2[i] = y[i];
188  }
189  else if( strcmp(plotCom.chPType[np],"OUTW") == 0 )
190  {
191  /* plot outward and attenuated incident continuum */
192  chSymPlt1 = 'i';
193  y[i] = (realnum)log10(MAX2(rfield.flux[i]*opac.opacity_abs[i],
194  1e-37));
195  strcpy( chSymPlt2, "o " );
196  y2[i] = (realnum)log10(MAX2((rfield.outlin[i]+rfield.outlin_noplot[i]+rfield.ConEmitOut[i])*
197  opac.opacity_abs[i],1e-37));
198  }
199  else if( strcmp(plotCom.chPType[np],"DIFF") == 0 )
200  {
201  /* plot "diffuse" continuum from last zone only */
202  y[i] = (realnum)log10(MAX2(rfield.ConEmitLocal[nzone][i]*rfield.anu2[i]/
203  rfield.widflx[i],1e-37));
204  y2[i] = y[i];
205  }
206  else if( strcmp(plotCom.chPType[np],"CRAW") == 0 )
207  {
208  y[i] = (realnum)log10(MAX2(rfield.flux_total_incident[i],1e-37));
209  y2[i] = (realnum)MAX2((rfield.flux[i]+
211  rfield.ConEmitOut[i])*radius.r1r0sq,1e-37);
212  y2[i] = (realnum)log10(y2[i]);
213  }
214 
215  if( y[i] > -36.9 )
216  {
217  ymin = MIN2(ymin,(double)y[i]);
218  }
219 
220  if( y2[i] > -36.9 )
221  {
222  ymin2 = MIN2(ymin2,(double)y2[i]);
223  }
224 
225  ymax = MAX2(ymax,(double)y[i]);
226  ymax = MAX2(ymax,(double)y2[i]);
227  }
228  }
229 
230  if( trace.lgTrace )
231  {
232  fprintf( ioQQQ, " PLOT called for the first time, YMAX, MIN=%10.2e%10.2e\n",
233  ymax, ymin );
234  }
235 
236  /* lower min by at most 5 dex below peak */
237  ymin2 = MAX3(ymax-5.,-35.,ymin2);
238 
239  /* make sure there is room at the bottom */
240  ymin = MIN3(ymin2,ymin,ymax-1.);
241 
242  /* emitted continuum is thermal, so goes to zero */
243  if( strcmp(plotCom.chPType[np],"EMIT") == 0 )
244  {
245  ymin = MAX2(ymin,ymax-4.);
246  }
247 
248  if( plotCom.lgPltTrace[np] )
249  {
250  fprintf( ioQQQ, " YMAX, MIN=%14.4e%14.4e Npnts=%4ld\n",
251  ymax
252  , ymin, rfield.nflux );
253  }
254  strcpy( chXtitle, "Log(nu fnu) vs LOG(nu)" );
255 
256  chSym = chSymPlt1;
257 
258  pltr(rfield.anulog,y,rfield.nflux,xmin,xmax,ymin,ymax,chSym,chXtitle
259  ,1,plotCom.lgPltTrace[np]);
260 
261  chSym = chSymPlt2[0];
262 
263  pltr(rfield.anulog,y2,rfield.nflux,xmin,xmax,ymin,ymax,chSym,chXtitle
264  ,3,plotCom.lgPltTrace[np]);
265 
266  free( y );
267  free( y2 );
268  return;
269 }
270 
271 /*pltmap generate plot of heating and cooling map */
272 
274  long int np,
275  const char *chCall)
276 {
277  char chXtitle[23];
278  static bool lgTlkSav;
279  char chSym;
280 
281  long int i;
282 
283  static double xmax,
284  xmin,
285  ymax,
286  ymin;
287 
288  DEBUG_ENTRY( "pltmap()" );
289 
290  if( strcmp(chCall,"FIRST") == 0 )
291  {
292  return;
293  }
294 
295  lgTlkSav = called.lgTalk;
296  called.lgTalk = false;
297  hcmap.lgMapBeingDone = true;
298  hcmap.RangeMap[0] = (realnum)pow((realnum)10.f,plotCom.pltxmn[np]);
299  hcmap.RangeMap[1] = (realnum)pow((realnum)10.f,plotCom.pltxmx[np]);
300  map_do(ioQQQ, " map");
301  called.lgTalk = lgTlkSav;
302 
303  for( i=0; i < hcmap.nmap; i++ )
304  {
305  hcmap.temap[i] = (realnum)log10(hcmap.temap[i]);
306  }
307 
308  xmin = MIN2(hcmap.temap[0],hcmap.temap[hcmap.nmap-1]);
309  xmin = MAX2((double)plotCom.pltxmn[np],xmin);
310  xmax = MAX2(hcmap.temap[0],hcmap.temap[hcmap.nmap-1]);
311  xmax = MIN2(xmax,(double)plotCom.pltxmx[np]);
312 
313  if( plotCom.lgPltTrace[np] )
314  {
315  fprintf( ioQQQ, " xmin, xmax=%12.4e%12.4e nmap=%4ld\n",
316  xmin, xmax, hcmap.nmap );
317  }
318 
319  ymin = FLT_MAX;
320  ymax = -FLT_MAX;
321 
322  for( i=0; i < hcmap.nmap; i++ )
323  {
324  if( (double)hcmap.temap[i] > xmin && (double)hcmap.temap[i] < xmax )
325  {
326  hcmap.hmap[i] = (realnum)log10(MAX2(hcmap.hmap[i],1e-35));
327  hcmap.cmap[i] = (realnum)log10(MAX2(hcmap.cmap[i],1e-35));
328  if( hcmap.cmap[i] > -34. )
329  {
330  ymin = MIN3(ymin,hcmap.hmap[i],hcmap.cmap[i]);
331  }
332  else
333  {
334  ymin = MIN2(ymin,(double)hcmap.hmap[i]);
335  }
336  ymax = MAX3(ymax,hcmap.hmap[i],hcmap.cmap[i]);
337  }
338  }
339 
340  if( trace.lgTrace )
341  {
342  fprintf( ioQQQ, " PLOT called for the first time, YMAX, MIN=%10.2e%10.2e\n",
343  ymax, ymin );
344  }
345 
346  if( plotCom.lgPltTrace[np] )
347  {
348  fprintf( ioQQQ, " YMAX, MIN=%14.4e%14.4e Npnts=%4ld\n",
349  ymax, ymin, hcmap.nmap );
350  }
351 
352  chSym = 'H';
353  strcpy( chXtitle, "heating - cooling v te" );
354 
355  pltr(hcmap.temap,hcmap.hmap,hcmap.nmap,xmin,xmax,ymin,ymax,chSym,
356  chXtitle,1,plotCom.lgPltTrace[np]);
357 
358  chSym = 'C';
359 
360  pltr(hcmap.temap,hcmap.cmap,hcmap.nmap,xmin,xmax,ymin,ymax,chSym,
361  chXtitle,3,plotCom.lgPltTrace[np]);
362 
363  return;
364 }
365 
366 /*pltopc generate plot of local gas opacity */
368  long int np,
369  const char *chCall)
370 {
371  char chXtitle[23];
372  char chSym;
373  long int i;
374  double arg1,
375  arg2;
376  static double xmax,
377  xmin,
378  ymax,
379  ymin;
380  static realnum *y/*[rfield.nupper]*/,
381  *y2/*[rfield.nupper]*/;
382 
383  DEBUG_ENTRY( "pltopc()" );
384 
385  if( strcmp(chCall,"FIRST") == 0 )
386  {
387  return;
388  }
389 
390  y = (realnum*)MALLOC((size_t)rfield.nupper*sizeof(realnum) );
391  y2 = (realnum*)MALLOC((size_t)rfield.nupper*sizeof(realnum) );
392 
393  xmin = rfield.anulog[0];
394  xmin = MAX2((double)plotCom.pltxmn[np],xmin);
395  xmax = rfield.anulog[rfield.nflux-1];
396  xmax = MIN2(xmax,(double)plotCom.pltxmx[np]);
397 
398  if( plotCom.lgPltTrace[np] )
399  {
400  fprintf( ioQQQ, " XMIN, XMAX=%12.4e%12.4e NFLUX=%4ld\n",
401  xmin, xmax, rfield.nflux );
402  }
403 
404  ymin = FLT_MAX;
405  ymax = -FLT_MAX;
406 
407  for( i=0; i < rfield.nflux; i++ )
408  {
409  if( strcmp(plotCom.chPType[np],"OPAA") == 0 )
410  {
411  /* absorption opacity */
412  arg1 = opac.opacity_abs_savzon1[i];
413  arg2 = opac.opacity_abs[i];
414  }
415 
416  else if( strcmp(plotCom.chPType[np],"OPAS") == 0 )
417  {
418  /* scattering opacity */
419  arg1 = opac.opacity_sct_savzon1[i];
420  arg2 = opac.opacity_sct[i];
421  }
422 
423  else if( strcmp(plotCom.chPType[np],"OPAT") == 0 )
424  {
425  /* total opacity */
427  arg2 = opac.opacity_abs[i] + opac.opacity_sct[i];
428  }
429 
430  else
431  {
432  /* this cannot happen since type was set to one of above */
433  fprintf( ioQQQ, " pltopc type=%4.4s not known. STOP\n",
434  plotCom.chPType[np] );
435  cdEXIT(EXIT_FAILURE);
436  }
437 
438  y[i] = (realnum)log10(MAX2(arg1/dense.gas_phase[ipHYDROGEN],1e-35));
439  y2[i] = (realnum)log10(MAX2(arg2/dense.gas_phase[ipHYDROGEN],1e-35));
440 
441  if( (double)rfield.anulog[i] > xmin && (double)rfield.anulog[i] < xmax )
442  {
443  ymin = MIN3(ymin,y[i],y2[i]);
444  ymax = MAX3(ymax,y[i],y2[i]);
445  }
446  }
447 
448  if( trace.lgTrace )
449  {
450  fprintf( ioQQQ, " PLOT called for the first time, YMAX, MIN=%10.2e%10.2e\n",
451  ymax, ymin );
452  }
453 
454  /* lower min by factor of 10 to show absorption in next plot */
455  ymin = MAX2(ymin-1.,-35.);
456  ymax += 1.;
457  if( plotCom.lgPltTrace[np] )
458  {
459  fprintf( ioQQQ, " YMAX, MIN=%14.4e%14.4e Npnts=%4ld\n",
460  ymax, ymin, rfield.nflux );
461  }
462 
463  strcpy( chXtitle, "Log(opacity) vs log(n)" );
464 
465  chSym = '.';
466  pltr(rfield.anulog,y,rfield.nflux,xmin,xmax,ymin,ymax,chSym,chXtitle
467  ,1,plotCom.lgPltTrace[np]);
468 
469  chSym = 'o';
470  pltr(rfield.anulog,y2,rfield.nflux,xmin,xmax,ymin,ymax,chSym,chXtitle
471  ,3,plotCom.lgPltTrace[np]);
472 
473  free(y);
474  free(y2);
475  return;
476 }
477 
478 /*pltr core plotting routine for generating line printer plots */
479 
481  /* the x-axis */
482  realnum x[],
483  /* the y-axi */
484  realnum y[],
485  /* number of points */
486  long int npnts,
487  /* mins and maxs, log of min and max of x-axis */
488  double xmin,
489  double xmax,
490  double ymin,
491  double ymax,
492  /* plot symbol */
493  char chSym,
494  char *chXtitle,
495  long int itim,
496  bool lgTrace)
497 {
498  static char chPage[59][122];
499 
500  long int i,
501  ix,
502  iy,
503  j,
504  nc;
505 
506  /* the max number of decades we can plot */
507 # define NDECAD 18
508 
509  static long int jpnt[NDECAD],
510  lowx,
511  lx;
512 
513  static double xdec,
514  xinc,
515  ydown,
516  yinc;
517 
518  /* this is log of smallestnumer in following set */
519  const realnum xAxisMin = -8.f;
520 
521  static char chLab[NDECAD][5]={"1E-8","1E-7","1E-6","1E-5",
522  "1E-4",".001","0.01"," 0.1"," 1 ",
523  " 10 "," 100","1000","1E4 ","1E5 ","1E6 ","1E7 ","1E8 ","1E9 "};
524 
525  DEBUG_ENTRY( "pltr()" );
526 
527  /* ITIM=1, first call, =2 intermediate calls, =3 for last call*/
528  if( itim == 1 )
529  {
530  /* first call, set left border of plot and clear out array */
531  for( i=1; i < IHI; i++ )
532  {
533  chPage[i][0] = 'l';
534  for( j=1; j < IWID; j++ )
535  {
536  chPage[i][j] = ' ';
537  }
538  }
539 
540  /* centered label for plot */
541  strcpy( chPage[1], " " );
542  strcat( chPage[1], chXtitle );
543  strcat( chPage[1], input.chTitle );
544 
545  /* one dex increments in x and y marked special */
546  i = 1;
547  ydown = 0.;
548  yinc = (realnum)(IHI-2)/(ymax - ymin);
549  nc = 0;
550 
551  while( i <= IHI && nc < 200 )
552  {
553  chPage[i-1][1] = '-';
554  ydown += yinc;
555  i = (long)(ydown + 1);
556  nc += 1;
557  }
558 
559  /* bottom increments of plot */
560  for( i=0; i < IWID; i++ )
561  {
562  chPage[IHI-1][i] = '-';
563  }
564 
565  if( xmin < xAxisMin )
566  {
567  fprintf(ioQQQ," plts: xmin is less than min value in array\n");
568  cdEXIT(EXIT_FAILURE);
569  }
570  /* LX is pointer to label for number in x-axis in chLab */
571  if( xmin < 0. )
572  {
573  lx = (long)(4.999-fabs(xmin));
574  /* lx is the offset within the array of x-axis values */
575  /* >>chng 99 jun 11 change to allow any min value of x-axis */
576  lx = (long)(fabs(xAxisMin)-0.001-fabs(xmin));
577  lx = MAX2(0,lx);
578  /* this is lowest decade on plot */
579  xdec = -floor(fabs(xmin)+1e-5);
580  }
581  else
582  {
583  double aa;
584  lx = (long)MAX2(0.,4.+xmin);
585  /* lx is the offset within the array of x-axis values */
586  lx = (long)MAX2(0.,4.+xmin);
587  /* >>chng 99 jun 11 change to allow any min value of x-axis */
588  aa = fabs(xAxisMin);
589  lx = (long)MAX2(0., aa-1. + xmin );
590  xdec = floor(xmin+1e-5);
591  }
592 
593  lowx = lx + 1;
594  xinc = (realnum)(IWID-1)/(xmax - xmin);
595  i = (long)MAX2(1.,(xdec-xmin)*xinc+1.);
596  nc = 0;
597 
598  while( i < IWID && nc < 100 )
599  {
600  chPage[IHI-2][i - 1] = 'l';
601 
602  /* fix position of x labels */
603  lx = MIN2(lx+1,NDECAD);
604 
605  /* slight offset to center label */
606  jpnt[lx-1] = MAX2(0,i-3);
607  jpnt[lx-1] = MIN2((long)IWID-4,jpnt[lx-1]);
608  xdec += 1.;
609  i = (long)MAX2(1.,(xdec-xmin)*xinc+1.);
610  nc += 1;
611  }
612  }
613 
614  /* everything falls down through here */
615  /* now fill in data, symbol is chSym */
616  for( i=0; i < npnts; i++ )
617  {
618  if( (double)x[i] > xmin && (double)x[i] < xmax )
619  {
620  iy = (long)(IHI - MAX2(y[i]-ymin,0.)*yinc);
621  iy = MAX2(1,iy);
622  ix = (long)((x[i] - xmin)*xinc + 1);
623 
624  if( lgTrace )
625  {
626  fprintf( ioQQQ, " x, y, ix, iy=%7.3f%7.3f%4ld%4ld\n",
627  x[i], y[i], ix, iy );
628  }
629  chPage[iy-1][ix - 1] = chSym;
630  }
631  }
632 
633  if( itim == 3 )
634  {
635  /* make the plot */
636  fprintf( ioQQQ, "1\n" );
637  for( i=1; i < IHI; i++ )
638  {
639  fprintf( ioQQQ, " %121.121s\n", chPage[i] );
640  }
641 
642  /* now put on label for X-axis */
643  for( i=0; i < IWID; i++ )
644  {
645  chPage[0][i] = ' ';
646  }
647 
648  for( i=lowx-1; i < lx; i++ )
649  {
650  /* copy the four char of the numeric string */
651  strncpy(chPage[0]+jpnt[i] , chLab[i+1] , 4);
652  }
653  fprintf( ioQQQ, " %121.121s\n", chPage[0] );
654  }
655  return;
656 }

Generated for cloudy by doxygen 1.8.3.1