00001
00002
00003
00004 #include "cddefines.h"
00005 #include "iterations.h"
00006 #include "called.h"
00007 #define IHI 59
00008 #define IWID 121
00009 #include "input.h"
00010 #include "rfield.h"
00011 #include "trace.h"
00012 #include "radius.h"
00013 #include "geometry.h"
00014 #include "opacity.h"
00015 #include "dense.h"
00016 #include "hcmap.h"
00017 #include "plot.h"
00018
00019 t_plotCom plotCom;
00020
00021
00022 STATIC void pltcon(long int np,
00023 const char *chCall);
00024
00025
00026 STATIC void pltmap(long int np,
00027 const char *chCall);
00028
00029
00030 STATIC void pltopc(long int np,
00031 const char *chCall);
00032
00033
00034 STATIC void pltr(realnum[],realnum[],long,double,double,double,double,
00035 char,char*,long,bool);
00036
00037
00038 void plot(const char *chCall)
00039 {
00040 long int np;
00041
00042 DEBUG_ENTRY( "plot()" );
00043
00044
00045
00046 if( !plotCom.lgPlotON || !called.lgTalk )
00047 {
00048 return;
00049 }
00050
00051 if( !iterations.lgLastIt && (strcmp(chCall,"FIRST") != 0) )
00052 {
00053 return;
00054 }
00055
00056
00057 for( np=0; np < plotCom.nplot; np++ )
00058 {
00059
00060 if( strcmp(plotCom.chPType[np]," MAP") == 0 )
00061 {
00062
00063 pltmap(np,chCall);
00064 }
00065 else if( strcmp(plotCom.chPType[np] ,"CONT") == 0 ||
00066 strcmp(plotCom.chPType[np] ,"CRAW") == 0 ||
00067 strcmp(plotCom.chPType[np] ,"DIFF") == 0 ||
00068 strcmp(plotCom.chPType[np] ,"REFL") == 0 ||
00069 strcmp(plotCom.chPType[np] ,"EMIT") == 0 ||
00070 strcmp(plotCom.chPType[np] ,"CPHT") == 0 ||
00071 strcmp(plotCom.chPType[np] ,"OUTW") == 0 )
00072 {
00073
00074 pltcon(np,chCall);
00075 }
00076
00077 else if(
00078 strcmp(plotCom.chPType[np] ,"OPAA") == 0 ||
00079 strcmp(plotCom.chPType[np] ,"OPAS") == 0 ||
00080 strcmp(plotCom.chPType[np] ,"OPAT") == 0 )
00081 {
00082
00083 pltopc(np,chCall);
00084 }
00085 else
00086 {
00087 fprintf( ioQQQ, " PLOT type=%4.4s not known. STOP\n",
00088 plotCom.chPType[np] );
00089 cdEXIT(EXIT_FAILURE);
00090 }
00091 }
00092
00093 return;
00094 }
00095
00096
00097
00098 STATIC void pltcon(
00099 long int np,
00100 const char *chCall)
00101 {
00102 char chSymPlt2[3],
00103 chXtitle[23];
00104 char chSym,
00105 chSymPlt1;
00106 long int i;
00107 double contin,
00108 ymin2;
00109 static double xmax,
00110 xmin,
00111 ymax,
00112 ymin;
00113 static realnum *y,
00114 *y2;
00115
00116 DEBUG_ENTRY( "pltcon()" );
00117
00118 if( strcmp(chCall,"FIRST") == 0 )
00119 {
00120 return;
00121 }
00122
00123 xmin = rfield.anulog[0];
00124 xmin = MAX2((double)plotCom.pltxmn[np],xmin);
00125 xmax = rfield.anulog[rfield.nflux-1];
00126 xmax = MIN2(xmax,(double)plotCom.pltxmx[np]);
00127
00128 if( plotCom.lgPltTrace[np] )
00129 {
00130 fprintf( ioQQQ, " XMIN, XMAX=%12.4e%12.4e NFLUX=%4ld\n",
00131 xmin, xmax, rfield.nflux );
00132 }
00133
00134 if( strcmp(plotCom.chPType[np],"REFL") == 0 && geometry.lgSphere )
00135 {
00136 fprintf( ioQQQ, " Reflected continuum not computed when SPHERE set.\n" );
00137 return;
00138 }
00139
00140 y = (realnum*)MALLOC((size_t)rfield.nupper*sizeof(realnum) );
00141 y2 = (realnum*)MALLOC((size_t)rfield.nupper*sizeof(realnum) );
00142
00143
00144 chSymPlt1 = '.';
00145 strcpy( chSymPlt2, "o " );
00146 ymin = FLT_MAX;
00147 ymin2 = FLT_MAX;
00148 ymax = -FLT_MAX;
00149 for( i=0; i < rfield.nflux; i++ )
00150 {
00151 if( (double)rfield.anulog[i] > xmin && (double)rfield.anulog[i] < xmax )
00152 {
00153 if( strcmp(plotCom.chPType[np],"CONT") == 0 )
00154 {
00155 y[i] = (realnum)log10(MAX2(rfield.flux_total_incident[0][i]/rfield.widflx[i]*
00156 rfield.anu2[i],1e-37));
00157
00158 contin = rfield.flux[0][i] + rfield.ConEmitOut[0][i]*geometry.covgeo + rfield.ConEmitReflec[0][i];
00159 y2[i] = (realnum)MAX2((contin/rfield.widflx[i]+(rfield.outlin[0][i]+rfield.outlin_noplot[i])/
00160 rfield.anu[i]*geometry.covgeo)*
00161 rfield.anu2[i]*radius.r1r0sq,1e-37);
00162 y2[i] = (realnum)log10(y2[i]);
00163 }
00164 else if( strcmp(plotCom.chPType[np],"CPHT") == 0 )
00165 {
00166
00167 y[i] = (realnum)log10(MAX2(rfield.flux_total_incident[0][i]/rfield.widflx[i],
00168 1e-37));
00169 contin = rfield.flux[0][i] + rfield.ConEmitOut[0][i]*geometry.covgeo/rfield.widflx[i];
00170 y2[i] = (realnum)MAX2((contin+(rfield.outlin[0][i]+rfield.outlin[0][i])/rfield.anu[i]*
00171 geometry.covgeo)* radius.r1r0sq,1e-37);
00172 y2[i] = (realnum)log10(y2[i]);
00173 }
00174 else if( strcmp(plotCom.chPType[np],"REFL") == 0 )
00175 {
00176
00177 y[i] = (realnum)log10(MAX2((rfield.ConEmitReflec[0][i]/rfield.widflx[i]+
00178 rfield.reflin[0][i])*rfield.anu2[i],1e-37));
00179 y2[i] = y[i];
00180 }
00181 else if( strcmp(plotCom.chPType[np],"EMIT") == 0 )
00182 {
00183
00184 y[i] = (realnum)log10(MAX2(
00185 ((rfield.ConEmitReflec[0][i]+rfield.ConEmitOut[0][i])/
00186 rfield.widflx[i]+
00187 (rfield.outlin[0][i]+rfield.outlin_noplot[i]+rfield.reflin[0][i])/rfield.anu[i] )*
00188 rfield.anu2[i],1e-37));
00189 y2[i] = y[i];
00190 }
00191 else if( strcmp(plotCom.chPType[np],"OUTW") == 0 )
00192 {
00193
00194 chSymPlt1 = 'i';
00195 y[i] = (realnum)log10(MAX2(rfield.flux[0][i]*opac.opacity_abs[i],
00196 1e-37));
00197 strcpy( chSymPlt2, "o " );
00198 y2[i] = (realnum)log10(MAX2((rfield.outlin[0][i]+rfield.outlin_noplot[i]+rfield.ConEmitOut[0][i])*
00199 opac.opacity_abs[i],1e-37));
00200 }
00201 else if( strcmp(plotCom.chPType[np],"DIFF") == 0 )
00202 {
00203
00204 y[i] = (realnum)log10(MAX2(rfield.ConEmitLocal[nzone][i]*rfield.anu2[i]/
00205 rfield.widflx[i],1e-37));
00206 y2[i] = y[i];
00207 }
00208 else if( strcmp(plotCom.chPType[np],"CRAW") == 0 )
00209 {
00210 y[i] = (realnum)log10(MAX2(rfield.flux_total_incident[0][i],1e-37));
00211 y2[i] = (realnum)MAX2((rfield.flux[0][i]+
00212 rfield.otscon[i]+rfield.otslin[i]+rfield.outlin[0][i]+rfield.outlin_noplot[i]+
00213 rfield.ConEmitOut[0][i])*radius.r1r0sq,1e-37);
00214 y2[i] = (realnum)log10(y2[i]);
00215 }
00216
00217 if( y[i] > -36.9 )
00218 {
00219 ymin = MIN2(ymin,(double)y[i]);
00220 }
00221
00222 if( y2[i] > -36.9 )
00223 {
00224 ymin2 = MIN2(ymin2,(double)y2[i]);
00225 }
00226
00227 ymax = MAX2(ymax,(double)y[i]);
00228 ymax = MAX2(ymax,(double)y2[i]);
00229 }
00230 }
00231
00232 if( trace.lgTrace )
00233 {
00234 fprintf( ioQQQ, " PLOT called for the first time, YMAX, MIN=%10.2e%10.2e\n",
00235 ymax, ymin );
00236 }
00237
00238
00239 ymin2 = MAX3(ymax-5.,-35.,ymin2);
00240
00241
00242 ymin = MIN3(ymin2,ymin,ymax-1.);
00243
00244
00245 if( strcmp(plotCom.chPType[np],"EMIT") == 0 )
00246 {
00247 ymin = MAX2(ymin,ymax-4.);
00248 }
00249
00250 if( plotCom.lgPltTrace[np] )
00251 {
00252 fprintf( ioQQQ, " YMAX, MIN=%14.4e%14.4e Npnts=%4ld\n",
00253 ymax
00254 , ymin, rfield.nflux );
00255 }
00256 strcpy( chXtitle, "Log(nu fnu) vs LOG(nu)" );
00257
00258 chSym = chSymPlt1;
00259
00260 pltr(rfield.anulog,y,rfield.nflux,xmin,xmax,ymin,ymax,chSym,chXtitle
00261 ,1,plotCom.lgPltTrace[np]);
00262
00263 chSym = chSymPlt2[0];
00264
00265 pltr(rfield.anulog,y2,rfield.nflux,xmin,xmax,ymin,ymax,chSym,chXtitle
00266 ,3,plotCom.lgPltTrace[np]);
00267
00268 free( y );
00269 free( y2 );
00270 return;
00271 }
00272
00273
00274
00275 STATIC void pltmap(
00276 long int np,
00277 const char *chCall)
00278 {
00279 char chXtitle[23];
00280 static bool lgTlkSav;
00281 char chSym;
00282
00283 long int i;
00284
00285 static double xmax,
00286 xmin,
00287 ymax,
00288 ymin;
00289
00290 DEBUG_ENTRY( "pltmap()" );
00291
00292 if( strcmp(chCall,"FIRST") == 0 )
00293 {
00294 return;
00295 }
00296
00297 lgTlkSav = called.lgTalk;
00298 called.lgTalk = false;
00299 hcmap.lgMapBeingDone = true;
00300 hcmap.RangeMap[0] = (realnum)pow((realnum)10.f,plotCom.pltxmn[np]);
00301 hcmap.RangeMap[1] = (realnum)pow((realnum)10.f,plotCom.pltxmx[np]);
00302 map_do(ioQQQ, " map");
00303 called.lgTalk = lgTlkSav;
00304
00305 for( i=0; i < hcmap.nmap; i++ )
00306 {
00307 hcmap.temap[i] = (realnum)log10(hcmap.temap[i]);
00308 }
00309
00310 xmin = MIN2(hcmap.temap[0],hcmap.temap[hcmap.nmap-1]);
00311 xmin = MAX2((double)plotCom.pltxmn[np],xmin);
00312 xmax = MAX2(hcmap.temap[0],hcmap.temap[hcmap.nmap-1]);
00313 xmax = MIN2(xmax,(double)plotCom.pltxmx[np]);
00314
00315 if( plotCom.lgPltTrace[np] )
00316 {
00317 fprintf( ioQQQ, " xmin, xmax=%12.4e%12.4e nmap=%4ld\n",
00318 xmin, xmax, hcmap.nmap );
00319 }
00320
00321 ymin = FLT_MAX;
00322 ymax = -FLT_MAX;
00323
00324 for( i=0; i < hcmap.nmap; i++ )
00325 {
00326 if( (double)hcmap.temap[i] > xmin && (double)hcmap.temap[i] < xmax )
00327 {
00328 hcmap.hmap[i] = (realnum)log10(MAX2(hcmap.hmap[i],1e-35));
00329 hcmap.cmap[i] = (realnum)log10(MAX2(hcmap.cmap[i],1e-35));
00330 if( hcmap.cmap[i] > -34. )
00331 {
00332 ymin = MIN3(ymin,hcmap.hmap[i],hcmap.cmap[i]);
00333 }
00334 else
00335 {
00336 ymin = MIN2(ymin,(double)hcmap.hmap[i]);
00337 }
00338 ymax = MAX3(ymax,hcmap.hmap[i],hcmap.cmap[i]);
00339 }
00340 }
00341
00342 if( trace.lgTrace )
00343 {
00344 fprintf( ioQQQ, " PLOT called for the first time, YMAX, MIN=%10.2e%10.2e\n",
00345 ymax, ymin );
00346 }
00347
00348 if( plotCom.lgPltTrace[np] )
00349 {
00350 fprintf( ioQQQ, " YMAX, MIN=%14.4e%14.4e Npnts=%4ld\n",
00351 ymax, ymin, hcmap.nmap );
00352 }
00353
00354 chSym = 'H';
00355 strcpy( chXtitle, "heating - cooling v te" );
00356
00357 pltr(hcmap.temap,hcmap.hmap,hcmap.nmap,xmin,xmax,ymin,ymax,chSym,
00358 chXtitle,1,plotCom.lgPltTrace[np]);
00359
00360 chSym = 'C';
00361
00362 pltr(hcmap.temap,hcmap.cmap,hcmap.nmap,xmin,xmax,ymin,ymax,chSym,
00363 chXtitle,3,plotCom.lgPltTrace[np]);
00364
00365 return;
00366 }
00367
00368
00369 STATIC void pltopc(
00370 long int np,
00371 const char *chCall)
00372 {
00373 char chXtitle[23];
00374 char chSym;
00375 long int i;
00376 double arg1,
00377 arg2;
00378 static double xmax,
00379 xmin,
00380 ymax,
00381 ymin;
00382 static realnum *y,
00383 *y2;
00384
00385 DEBUG_ENTRY( "pltopc()" );
00386
00387 if( strcmp(chCall,"FIRST") == 0 )
00388 {
00389 return;
00390 }
00391
00392 y = (realnum*)MALLOC((size_t)rfield.nupper*sizeof(realnum) );
00393 y2 = (realnum*)MALLOC((size_t)rfield.nupper*sizeof(realnum) );
00394
00395 xmin = rfield.anulog[0];
00396 xmin = MAX2((double)plotCom.pltxmn[np],xmin);
00397 xmax = rfield.anulog[rfield.nflux-1];
00398 xmax = MIN2(xmax,(double)plotCom.pltxmx[np]);
00399
00400 if( plotCom.lgPltTrace[np] )
00401 {
00402 fprintf( ioQQQ, " XMIN, XMAX=%12.4e%12.4e NFLUX=%4ld\n",
00403 xmin, xmax, rfield.nflux );
00404 }
00405
00406 ymin = FLT_MAX;
00407 ymax = -FLT_MAX;
00408
00409 for( i=0; i < rfield.nflux; i++ )
00410 {
00411 if( strcmp(plotCom.chPType[np],"OPAA") == 0 )
00412 {
00413
00414 arg1 = opac.opacity_abs_savzon1[i];
00415 arg2 = opac.opacity_abs[i];
00416 }
00417
00418 else if( strcmp(plotCom.chPType[np],"OPAS") == 0 )
00419 {
00420
00421 arg1 = opac.opacity_sct_savzon1[i];
00422 arg2 = opac.opacity_sct[i];
00423 }
00424
00425 else if( strcmp(plotCom.chPType[np],"OPAT") == 0 )
00426 {
00427
00428 arg1 = opac.opacity_abs_savzon1[i] + opac.opacity_sct_savzon1[i];
00429 arg2 = opac.opacity_abs[i] + opac.opacity_sct[i];
00430 }
00431
00432 else
00433 {
00434
00435 fprintf( ioQQQ, " pltopc type=%4.4s not known. STOP\n",
00436 plotCom.chPType[np] );
00437 cdEXIT(EXIT_FAILURE);
00438 }
00439
00440 y[i] = (realnum)log10(MAX2(arg1/dense.gas_phase[ipHYDROGEN],1e-35));
00441 y2[i] = (realnum)log10(MAX2(arg2/dense.gas_phase[ipHYDROGEN],1e-35));
00442
00443 if( (double)rfield.anulog[i] > xmin && (double)rfield.anulog[i] < xmax )
00444 {
00445 ymin = MIN3(ymin,y[i],y2[i]);
00446 ymax = MAX3(ymax,y[i],y2[i]);
00447 }
00448 }
00449
00450 if( trace.lgTrace )
00451 {
00452 fprintf( ioQQQ, " PLOT called for the first time, YMAX, MIN=%10.2e%10.2e\n",
00453 ymax, ymin );
00454 }
00455
00456
00457 ymin = MAX2(ymin-1.,-35.);
00458 ymax += 1.;
00459 if( plotCom.lgPltTrace[np] )
00460 {
00461 fprintf( ioQQQ, " YMAX, MIN=%14.4e%14.4e Npnts=%4ld\n",
00462 ymax, ymin, rfield.nflux );
00463 }
00464
00465 strcpy( chXtitle, "Log(opacity) vs log(n)" );
00466
00467 chSym = '.';
00468 pltr(rfield.anulog,y,rfield.nflux,xmin,xmax,ymin,ymax,chSym,chXtitle
00469 ,1,plotCom.lgPltTrace[np]);
00470
00471 chSym = 'o';
00472 pltr(rfield.anulog,y2,rfield.nflux,xmin,xmax,ymin,ymax,chSym,chXtitle
00473 ,3,plotCom.lgPltTrace[np]);
00474
00475 free(y);
00476 free(y2);
00477 return;
00478 }
00479
00480
00481
00482 STATIC void pltr(
00483
00484 realnum x[],
00485
00486 realnum y[],
00487
00488 long int npnts,
00489
00490 double xmin,
00491 double xmax,
00492 double ymin,
00493 double ymax,
00494
00495 char chSym,
00496 char *chXtitle,
00497 long int itim,
00498 bool lgTrace)
00499 {
00500 static char chPage[59][122];
00501
00502 long int i,
00503 ix,
00504 iy,
00505 j,
00506 nc;
00507
00508
00509 # define NDECAD 18
00510
00511 static long int jpnt[NDECAD],
00512 lowx,
00513 lx;
00514
00515 static double xdec,
00516 xinc,
00517 ydown,
00518 yinc;
00519
00520
00521 const realnum xAxisMin = -8.f;
00522
00523 static char chLab[NDECAD][5]={"1E-8","1E-7","1E-6","1E-5",
00524 "1E-4",".001","0.01"," 0.1"," 1 ",
00525 " 10 "," 100","1000","1E4 ","1E5 ","1E6 ","1E7 ","1E8 ","1E9 "};
00526
00527 DEBUG_ENTRY( "pltr()" );
00528
00529
00530 if( itim == 1 )
00531 {
00532
00533 for( i=1; i < IHI; i++ )
00534 {
00535 chPage[i][0] = 'l';
00536 for( j=1; j < IWID; j++ )
00537 {
00538 chPage[i][j] = ' ';
00539 }
00540 }
00541
00542
00543 strcpy( chPage[1], " " );
00544 strcat( chPage[1], chXtitle );
00545 strcat( chPage[1], input.chTitle );
00546
00547
00548 i = 1;
00549 ydown = 0.;
00550 yinc = (realnum)(IHI-2)/(ymax - ymin);
00551 nc = 0;
00552
00553 while( i <= IHI && nc < 200 )
00554 {
00555 chPage[i-1][1] = '-';
00556 ydown += yinc;
00557 i = (long)(ydown + 1);
00558 nc += 1;
00559 }
00560
00561
00562 for( i=0; i < IWID; i++ )
00563 {
00564 chPage[IHI-1][i] = '-';
00565 }
00566
00567 if( xmin < xAxisMin )
00568 {
00569 fprintf(ioQQQ," plts: xmin is less than min value in array\n");
00570 cdEXIT(EXIT_FAILURE);
00571 }
00572
00573 if( xmin < 0. )
00574 {
00575 lx = (long)(4.999-fabs(xmin));
00576
00577
00578 lx = (long)(fabs(xAxisMin)-0.001-fabs(xmin));
00579 lx = MAX2(0,lx);
00580
00581 xdec = -floor(fabs(xmin)+1e-5);
00582 }
00583 else
00584 {
00585 double aa;
00586 lx = (long)MAX2(0.,4.+xmin);
00587
00588 lx = (long)MAX2(0.,4.+xmin);
00589
00590 aa = fabs(xAxisMin);
00591 lx = (long)MAX2(0., aa-1. + xmin );
00592 xdec = floor(xmin+1e-5);
00593 }
00594
00595 lowx = lx + 1;
00596 xinc = (realnum)(IWID-1)/(xmax - xmin);
00597 i = (long)MAX2(1.,(xdec-xmin)*xinc+1.);
00598 nc = 0;
00599
00600 while( i < IWID && nc < 100 )
00601 {
00602 chPage[IHI-2][i - 1] = 'l';
00603
00604
00605 lx = MIN2(lx+1,NDECAD);
00606
00607
00608 jpnt[lx-1] = MAX2(0,i-3);
00609 jpnt[lx-1] = MIN2((long)IWID-4,jpnt[lx-1]);
00610 xdec += 1.;
00611 i = (long)MAX2(1.,(xdec-xmin)*xinc+1.);
00612 nc += 1;
00613 }
00614 }
00615
00616
00617
00618 for( i=0; i < npnts; i++ )
00619 {
00620 if( (double)x[i] > xmin && (double)x[i] < xmax )
00621 {
00622 iy = (long)(IHI - MAX2(y[i]-ymin,0.)*yinc);
00623 iy = MAX2(1,iy);
00624 ix = (long)((x[i] - xmin)*xinc + 1);
00625
00626 if( lgTrace )
00627 {
00628 fprintf( ioQQQ, " x, y, ix, iy=%7.3f%7.3f%4ld%4ld\n",
00629 x[i], y[i], ix, iy );
00630 }
00631 chPage[iy-1][ix - 1] = chSym;
00632 }
00633 }
00634
00635 if( itim == 3 )
00636 {
00637
00638 fprintf( ioQQQ, "1\n" );
00639 for( i=1; i < IHI; i++ )
00640 {
00641 fprintf( ioQQQ, " %121.121s\n", chPage[i] );
00642 }
00643
00644
00645 for( i=0; i < IWID; i++ )
00646 {
00647 chPage[0][i] = ' ';
00648 }
00649
00650 for( i=lowx-1; i < lx; i++ )
00651 {
00652
00653 strncpy(chPage[0]+jpnt[i] , chLab[i+1] , 4);
00654 }
00655 fprintf( ioQQQ, " %121.121s\n", chPage[0] );
00656 }
00657 return;
00658 }