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