cloudy: /home66/gary/public_html/cloudy/c08_branch/source/hcmap.cpp Source File

00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
00002  * others.  For conditions of distribution and use see copyright notice in license.txt */
00003 /*map_do produce map of heating-cooling space for specified zone, called as result of
00004  * map command  */
00005 #include "cddefines.h"
00006 #include "thermal.h"
00007 #include "cooling.h"
00008 #include "called.h"
00009 #include "dense.h"
00010 #include "mole.h"
00011 #include "phycon.h"
00012 #include "trace.h"
00013 #include "pressure.h"
00014 #include "conv.h"
00015 #include "hcmap.h"
00016 #ifdef EPS
00017 #       undef EPS
00018 #endif
00019 #define EPS     0.005
00020 
00021 void map_do(
00022                         FILE *io, 
00023                         const char *chType)
00024 {
00025         char chLabel[NCOLNT_LAB_LEN+1];
00026         char units;
00027         long int i, 
00028           ksav, 
00029           j, 
00030           jsav, 
00031           k,
00032           nelem;
00033         realnum wl;
00034         double cfrac, 
00035           ch, 
00036           fac, 
00037           factor, 
00038           hfrac, 
00039           oldch, 
00040           ratio, 
00041           strhet, 
00042           strong, 
00043           t1, 
00044           tlowst, 
00045           tmax, 
00046           tmin, 
00047           torg;
00048 
00049         DEBUG_ENTRY( "map_do()" );
00050 
00051         t1 = phycon.te;
00052         torg = phycon.te;
00053         hcmap.lgMapOK = true;
00054         /* flag indicating that we have computed a map */
00055         hcmap.lgMapDone = true;
00056 
00057         /* make sure pressure has been evaluated */
00058         /* this sets values of pressure.PresTotlCurr */
00059         PresTotCurrent();
00060 
00061         /* print out all coolants if all else fails */
00062         if( called.lgTalk )
00063         {
00064                 fprintf( io, " Cloudy punts, Te=%10.3e HTOT=%10.3e CTOT=%10.3e nzone=%4ld\n", 
00065                   phycon.te, thermal.htot, thermal.ctot, nzone );
00066                 fprintf( io, " COOLNG array is\n" );
00067 
00068                 if( thermal.ctot > 0. )
00069                 {
00070                         coolpr(io, "ZERO",1,0.,"ZERO");
00071                         for( i=0; i < thermal.ncltot; i++ )
00072                         {
00073                                 ratio = thermal.cooling[i]/thermal.ctot;
00074                                 if( ratio>EPS )
00075                                 {
00076                                         coolpr(io, (char*)thermal.chClntLab[i],thermal.collam[i],
00077                                           ratio,"DOIT");
00078                                 }
00079                         }
00080 
00081                         coolpr(io, "DONE",1,0.,"DONE");
00082                         fprintf( io, " Line heating array follows\n" );
00083                         coolpr(io, "ZERO",1,0.,"ZERO");
00084 
00085                         for( i=0; i < thermal.ncltot; i++ )
00086                         {
00087                                 ratio = thermal.heatnt[i]/thermal.ctot;
00088                                 if( ratio>EPS )
00089                                 {
00090                                         coolpr(io, (char*)thermal.chClntLab[i],thermal.collam[i],
00091                                           ratio,"DOIT");
00092                                 }
00093                         }
00094 
00095                         coolpr(io,"DONE",1,0.,"DONE");
00096                 }
00097         }
00098 
00099         /* map out te-ionization-cooling space before punching out. */
00100         if( called.lgTalk )
00101         {
00102                 fprintf( io, " map of heating, cooling, vs temp, follows.\n");
00103                 fprintf( io, 
00104                         "     Te     Heat--------------------> Cool--------------------->   dH/dT     dC/DT       Ne         NH      HII       Helium \n" );
00105         }
00106 
00107         if( strcmp(chType,"punt") == 0 )
00108         {
00109                 /* this is the original use of punt, we are punting
00110                  * only map within factor of two of final temperature
00111                  * fac will the range to either side of punted temperature */
00112                 fac = 1.5;
00113                 tmin = torg/fac;
00114                 tmax = torg*fac;
00115 
00116                 /* we want about 20 steps between high and low temperature
00117                  * default of nMapStep is 20, set with set nmaps command */
00118                 factor = pow(10.,log10(tmax/tmin)/(double)(hcmap.nMapStep));
00119                 TempChange(tmin , false);
00120         }
00121 
00122         else if( strcmp(chType," map") == 0 )
00123         {
00124                 /* create some sort of map of heating-cooling */
00125                 tlowst = MAX2(hcmap.RangeMap[0],phycon.TEMP_LIMIT_LOW);
00126                 tmin = tlowst*0.998;
00127                 tmax = MIN2(hcmap.RangeMap[1],phycon.TEMP_LIMIT_HIGH)*1.002;
00128 
00129                 /* we want about nMapStep (=20) steps between high and low temperature */
00130                 factor = pow(10.,log10(tmax/tmin)/(double)(hcmap.nMapStep));
00131                 double TeNew;
00132                 if( thermal.lgTeHigh )
00133                 {
00134                         /* high te */
00135                         factor = 1./factor;
00136                         /* TeHighest is highest possible temperature, 1E10 */
00137                         TeNew = (MIN2(hcmap.RangeMap[1],phycon.TEMP_LIMIT_HIGH)/factor);
00138                 }
00139 
00140                 else
00141                 {
00142                         /* low te */
00143                         TeNew = (tlowst/factor);
00144                 }
00145                 TempChange(TeNew , false);
00146         }
00147 
00148         else
00149         {
00150                 /* don't know what to do */
00151                 fprintf( ioQQQ, " PUNT called with insane argument,=%4.4s\n", 
00152                   chType );
00153                 cdEXIT(EXIT_FAILURE);
00154         }
00155 
00156         /* now allocate space for te, c, h vectors in map, if not already done */
00157         if( hcmap.nMapAlloc==0 )
00158         {
00159                 /* space not allocated, do so now */
00160                 hcmap.nMapAlloc = hcmap.nMapStep+4;
00161 
00162                 /* now make the space */
00163                 hcmap.temap = (realnum *)MALLOC( sizeof(realnum)*(size_t)(hcmap.nMapStep+4) );
00164                 if( hcmap.temap == NULL )
00165                 { 
00166                         printf( " not enough memory to allocate hcmap.temap in map_do\n" );
00167                         cdEXIT(EXIT_FAILURE);
00168                 }
00169                 hcmap.cmap = (realnum *)MALLOC( sizeof(realnum)*(size_t)(hcmap.nMapStep+4) );
00170                 if( hcmap.cmap == NULL )
00171                 { 
00172                         printf( " not enough memory to allocate hcmap.cmap in map_do\n" );
00173                         cdEXIT(EXIT_FAILURE);
00174                 }
00175                 hcmap.hmap = (realnum *)MALLOC( sizeof(realnum)*(size_t)(hcmap.nMapStep+4) );
00176                 if( hcmap.hmap == NULL )
00177                 { 
00178                         printf( " not enough memory to allocate hcmap.hmap in map_do\n" );
00179                         cdEXIT(EXIT_FAILURE);
00180                 }
00181         }
00182 
00183         thermal.lgCNegChk = false;
00184         hcmap.nmap = 0;
00185         oldch = 0.;
00186         TempChange(phycon.te *factor , true);
00187         if( trace.nTrConvg )
00188                 fprintf(ioQQQ, "    MAP called temp range %.4e %.4e in %li stops ===============================================\n",
00189                 tmin,
00190                 tmax,
00191                 hcmap.nmap);
00192 
00193         while( (double)phycon.te < tmax*0.999 && (double)phycon.te > tmin*1.001 )
00194         {
00195                 /* this sets values of pressure.PresTotlCurr */
00196                 PresTotCurrent();
00197 
00198                 /* must reset upper and lower bounds for ionization distributions */
00199                 /* fix number of stages of ionization */
00200                 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
00201                 {
00202                         if( dense.lgElmtOn[nelem] )
00203                         {
00204                                 dense.IonLow[nelem] = 0;
00205                                 /* 
00206                                 * IonHigh[n] is the highest stage of ionization present
00207                                 * the IonHigh array index is on the C scake, so [0] is hydrogen
00208                                 * the value is also on the C scale, so element [nelem] can range
00209                                 * from 0 to nelem+1 
00210                                 */
00211                                 dense.IonHigh[nelem] = nelem + 1;
00212                         }
00213                         else
00214                         {
00215                                 /* this element is turned off, make stages impossible */
00216                                 dense.IonLow[nelem] = -1;
00217                                 dense.IonHigh[nelem] = -1;
00218                         }
00219                 }
00220 
00221                 /* this turns on constant reevaluation of everything */
00222                 conv.lgSearch = true;
00223 
00224                 if( trace.nTrConvg )
00225                         fprintf(ioQQQ, "    MAP new temp %.4e ===============================================\n",
00226                         phycon.te );
00227 
00228                 /* this counts how many times ionize is called in this model after startr,
00229                  * and is flag used by ionize to understand it is being called the first time*/
00230                 conv.nTotalIoniz = 0;
00231 
00232                 /* this will force reinitialization of co network */
00233                 co.iteration_co = -2;
00234 
00235                 /* now get ionization solution for this temperature */
00236                 if( ConvEdenIoniz() )
00237                         lgAbort = true;
00238 
00239                 /* save results for later prints */
00240                 hcmap.temap[hcmap.nmap] = (realnum)phycon.te;
00241                 hcmap.cmap[hcmap.nmap] = (realnum)thermal.ctot;
00242                 hcmap.hmap[hcmap.nmap] = (realnum)thermal.htot;
00243 
00244                 wl = 0.f;
00245                 strong = 0.;
00246 
00247                 for( j=0; j < thermal.ncltot; j++ )
00248                 {
00249                         if( thermal.cooling[j] > strong )
00250                         {
00251                                 strcpy( chLabel, thermal.chClntLab[j] );
00252                                 strong = thermal.cooling[j];
00253                                 wl = thermal.collam[j];
00254                         }
00255                 }
00256 
00257                 cfrac = strong/thermal.ctot;
00258                 strhet = 0.;
00259                 /* these will be reset in following loop*/
00260                 ksav = -INT_MAX;
00261                 jsav = -INT_MAX;
00262 
00263                 for( k=0; k < LIMELM; k++ )
00264                 {
00265                         for( j=0; j < LIMELM; j++ )
00266                         {
00267                                 if( thermal.heating[k][j] > strhet )
00268                                 {
00269                                         strhet = thermal.heating[k][j];
00270                                         jsav = j;
00271                                         ksav = k;
00272                                 }
00273                         }
00274                 }
00275 
00276                 ch = thermal.ctot - thermal.htot;
00277                 /* use ratio to check for change of sign since product
00278                  * can underflow at low densities */
00279                 if( oldch/ch < 0. && called.lgTalk )
00280                 {
00281                         fprintf( io, " ----------------------------------------------- Probable thermal solution here. --------------------------------------------\n" );
00282                 }
00283 
00284                 oldch = ch;
00285                 hfrac = strhet/thermal.htot;
00286                 if( called.lgTalk )
00287                 {
00288                         /* convert to micros if IR transition */
00289                         if( wl < 100000.f )
00290                         {
00291                                 units = 'A';
00292                         }
00293 
00294                         else
00295                         {
00296                                 wl /= 10000.f;
00297                                 units = 'm';
00298                         }
00299 
00300                         if( trace.lgTrace )
00301                         {
00302                                 fprintf( io, "TRACE: te, htot, ctot%11.3e%11.3e%11.3e\n", 
00303                                   phycon.te, thermal.htot, thermal.ctot );
00304                         }
00305 
00306                         /*fprintf( io, 
00307                                 "%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",;*/
00308                         fprintf(io,"%s", PrintEfmt("%11.4e", phycon.te ) );
00309                         fprintf(io,"%s", PrintEfmt("%11.4e", thermal.htot ) );
00310                         fprintf(io," [%2ld][%2ld]%6.3f",
00311                           ksav, jsav,  
00312                           hfrac);
00313                         fprintf(io,"%s", PrintEfmt("%11.4e", thermal.ctot ) );
00314                         fprintf(io," %s %.1f%c%6.3f",
00315                           chLabel , 
00316                           wl, 
00317                           units, 
00318                           cfrac );
00319                         fprintf(io,"%s", PrintEfmt(" %10.2e", thermal.dHeatdT ) );
00320                         fprintf(io,"%s", PrintEfmt("%11.2e", thermal.dCooldT ) );
00321                         fprintf(io,"%s", PrintEfmt("%11.4e", dense.eden ) );
00322                         fprintf(io,"%s", PrintEfmt("%11.4e", dense.gas_phase[ipHYDROGEN] ) );
00323                         if( dense.lgElmtOn[ipHELIUM] )
00324                         {
00325                                 fprintf(io,"%6.2f%6.2f%6.2f%6.2f\n",
00326                                 log10(MAX2(1e-9,dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN])), 
00327                                 log10(MAX2(1e-9,dense.xIonDense[ipHELIUM][0]/dense.gas_phase[ipHELIUM])), 
00328                                 log10(MAX2(1e-9,dense.xIonDense[ipHELIUM][1]/dense.gas_phase[ipHELIUM])), 
00329                                 log10(MAX2(1e-9,dense.xIonDense[ipHELIUM][2]/dense.gas_phase[ipHELIUM])) );
00330                         }
00331                         fflush(io);
00332                 }
00333 
00334                 TempChange(phycon.te*factor , true);
00335                 /* increment nmap but do not exceed nMapAlloc */
00336                 hcmap.nmap = MIN2(hcmap.nMapAlloc,hcmap.nmap+1);
00337 
00338                 {
00339                         enum {DEBUG_LOC=false};
00340                         if( DEBUG_LOC )
00341                         {
00342                                 static int kount = 0;
00343                                 factor = 1.;
00344                                 TempChange(8674900. , true);
00345                                 ++kount;
00346                                 if( kount >=100 )
00347                                 {
00348                                         fprintf(ioQQQ," exiting in map_do\n");
00349                                         break;
00350                                 }
00351                         }
00352                 }
00353         }
00354 
00355         /* now check whether sharp inflections occurred, and also find the biggest jump
00356          * in the heating and cooling */
00357         hcmap.lgMapOK = true;
00358         /* >>chng 02 mar 04, lower bound had been 1, so [i-2] below was negative */
00359          for( i=2; i< hcmap.nmap-2; ++i )
00360          {
00361                 realnum s1,s2,s3;/* the three slopes we will use */
00362                 s1 = hcmap.cmap[i-2] - hcmap.cmap[i-1];
00363                 s2 = hcmap.cmap[i-1] - hcmap.cmap[i];
00364                 s3 = hcmap.cmap[i] - hcmap.cmap[i+1];
00365                 if( s1*s3 > 0. && s2*s3 < 0. )
00366                 {
00367                          /* of the three points, the outer had the same slope 
00368                           * (their product was positive) but there was an inflection
00369                           * between them (the negative product).  The data chain looked like
00370                           *     2 4
00371                           *    1 3  or vice versa, either case is wrong,
00372                           * with the logic in the above test, the problem point will aways be s2 */
00373                         fprintf( io,
00374                                 " cooling curve had double inflection at T=%.2e.  ",
00375                                 hcmap.temap[i]);
00376                         fprintf( io,    " Slopes were %.2e %.2e %.2e",  s1, s2, s3);
00377                         if( fabs(s2)/hcmap.cmap[i] > 0.05 )
00378                         {
00379                                 fprintf( io,
00380                                         " error large, (rel slope of %.2e).\n",
00381                                         s2 / hcmap.cmap[i]);
00382                                 hcmap.lgMapOK = false;
00383                         }
00384                         else
00385                         {
00386                                 fprintf( io,
00387                                         " error is small, (rel slope of %.2e).\n",
00388                                         s2 / hcmap.cmap[i]);
00389                         }
00390                 }
00391 
00392                 s1 = hcmap.hmap[i-2] - hcmap.hmap[i-1];
00393                 s2 = hcmap.hmap[i-1] - hcmap.hmap[i];
00394                 s3 = hcmap.hmap[i]   - hcmap.hmap[i+1];
00395                 if( s1*s3 > 0. && s2*s3 < 0. )
00396                 {
00397                          /* of the three points, the outer had the same slope 
00398                           * (their product was positive) but there was an inflection
00399                           * between them (the negative product).  The data chain looked like
00400                           *     2 4
00401                           *    1 3  or vice versa, either case is wrong */
00402                         fprintf( io,
00403                                 " heating curve had double inflection at T=%.2e.\n",
00404                                 hcmap.temap[i] );
00405                         hcmap.lgMapOK = false;
00406                 }
00407          }
00408 
00409         thermal.lgCNegChk = true;
00410         TempChange(t1 , false);
00411         return;
00412 }