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