00001
00002
00003
00004 #include "cddefines.h"
00005 #include "thermal.h"
00006 #include "radius.h"
00007 #include "conv.h"
00008 #include "lines_service.h"
00009 #include "dense.h"
00010 #include "taulines.h"
00011 #include "phycon.h"
00012 #include "elementnames.h"
00013 #include "dynamics.h"
00014 #include "punch.h"
00015
00016
00017
00018 static const int IPRINT= 9;
00019
00020
00021
00022 void HeatPunch(FILE* io)
00023 {
00024 char **chLabel,
00025 chLbl[11];
00026 bool lgHeatLine;
00027 int nFail;
00028 long int i,
00029 ipStrong,
00030 ipnt,
00031 *ipOrdered,
00032 *ipsave,
00033 j,
00034 *jpsave,
00035 k,
00036 level;
00037 double CS,
00038 ColHeat,
00039 EscP,
00040 Pump,
00041 Strong,
00042 TauIn,
00043 cool_total,
00044 heat_total;
00045 float *save;
00046
00047 DEBUG_ENTRY( "HeatPunch()" );
00048
00049 save = (float *) CALLOC(LIMELM*LIMELM, sizeof(float));
00050 ipsave = (long int *) CALLOC(LIMELM*LIMELM, sizeof(long int));
00051 jpsave = (long int *) CALLOC(LIMELM*LIMELM, sizeof(long int));
00052 ipOrdered = (long int *) CALLOC(LIMELM*LIMELM, sizeof(long int));
00053 chLabel = (char **) CALLOC(LIMELM*LIMELM, sizeof(char *));
00054
00055 for( i=0; i<LIMELM*LIMELM; ++i )
00056 {
00057 ipsave[i] = INT_MIN;
00058 jpsave[i] = INT_MIN;
00059 save[i] = -FLT_MAX;
00060 chLabel[i] = (char *) CALLOC(10, sizeof(char));
00061 }
00062
00063 cool_total = thermal.ctot;
00064 heat_total = thermal.htot;
00065
00066
00067
00068
00069
00070
00071 cool_total -= dynamics.Cool;
00072 heat_total -= dynamics.Heat;
00073 # if 0
00074 if(dynamics.Cool > dynamics.Heat)
00075 {
00076 cool_total -= dynamics.Heat;
00077 heat_total -= dynamics.Heat;
00078 }
00079 else
00080 {
00081 cool_total -= dynamics.Cool;
00082 heat_total -= dynamics.Cool;
00083 }
00084 # endif
00085
00086 ipnt = 0;
00087
00088
00089
00090
00091 for( i=0; i < LIMELM; i++ )
00092 {
00093 for( j=0; j < LIMELM; j++ )
00094 {
00095 if( thermal.heating[i][j]/heat_total > SMALLFLOAT )
00096 {
00097 ipsave[ipnt] = i;
00098 jpsave[ipnt] = j;
00099 save[ipnt] = (float)(thermal.heating[i][j]/heat_total);
00100 ipnt++;
00101 }
00102 }
00103 }
00104
00105
00106
00107
00108 for( i=0; i < thermal.ncltot; i++ )
00109 {
00110 if( thermal.heatnt[i]/heat_total > punch.WeakHeatCool )
00111 {
00112 float awl;
00113 awl = thermal.collam[i];
00114
00115 if( awl > 100000 )
00116 awl /= 10000;
00117 fprintf( io, " Negative coolant was %s %.2f %.2e\n",
00118 thermal.chClntLab[i], awl, thermal.heatnt[i]/heat_total );
00119 }
00120 }
00121
00122 if( !conv.lgConvTemp )
00123 {
00124 fprintf( io, "#>>>> Temperature not converged.\n" );
00125 }
00126 else if( !conv.lgConvEden )
00127 {
00128 fprintf( io, "#>>>> Electron density not converged.\n" );
00129 }
00130 else if( !conv.lgConvIoniz )
00131 {
00132 fprintf( io, "#>>>> Ionization not converged.\n" );
00133 }
00134 else if( !conv.lgConvPres )
00135 {
00136 fprintf( io, "#>>>> Pressure not converged.\n" );
00137 }
00138
00139
00140
00141 i = INT_MIN;
00142 j = INT_MIN;
00143
00144 for( k=0; k < ipnt; k++ )
00145 {
00146
00147
00148
00149 i = ipsave[k];
00150 j = jpsave[k];
00151
00152 if( i >= j )
00153 {
00154 if( dense.xIonDense[i][j] == 0. && thermal.heating[i][j]>SMALLFLOAT )
00155 fprintf(ioQQQ,"DISASTER assert about to be thrown - search for hit it\n");
00156
00157
00158
00159 ASSERT( dense.xIonDense[i][j] > 0. || thermal.heating[i][j]<SMALLFLOAT );
00160
00161 strcpy( chLabel[k], elementnames.chElementSym[i] );
00162 strcat( chLabel[k], elementnames.chIonStage[j] );
00163 }
00164
00165 else if( i == 0 && j == 1 )
00166 {
00167
00168
00169 strcpy( chLabel[k], "Hn=2" );
00170 }
00171 else if( i == 0 && j == 3 )
00172 {
00173
00174
00175 strcpy( chLabel[k], "Hion" );
00176 }
00177 else if( i == 0 && j == 7 )
00178 {
00179
00180 strcpy( chLabel[k], " UTA" );
00181 }
00182 else if( i == 0 && j == 8 )
00183 {
00184
00185
00186
00187 strcpy( chLabel[k], "H2vH" );
00188 }
00189 else if( i == 0 && j == 17 )
00190 {
00191
00192
00193
00194 strcpy( chLabel[k], "H2dH" );
00195 }
00196 else if( i == 0 && j == 9 )
00197 {
00198
00199 strcpy( chLabel[k], "COds" );
00200 }
00201 else if( i == 0 && j == 20 )
00202 {
00203
00204 strcpy( chLabel[k], "extH" );
00205 }
00206 else if( i == 0 && j == 11 )
00207 {
00208
00209 strcpy( chLabel[k], "H FF" );
00210 }
00211 else if( i == 0 && j == 12 )
00212 {
00213
00214 strcpy( chLabel[k], "Clin" );
00215 }
00216 else if( i == 0 && j == 13 )
00217 {
00218
00219 strcpy( chLabel[k], "GrnP" );
00220 }
00221 else if( i == 0 && j == 14 )
00222 {
00223
00224 strcpy( chLabel[k], "GrnC" );
00225 }
00226 else if( i == 0 && j == 15 )
00227 {
00228
00229 strcpy( chLabel[k], "H- " );
00230 }
00231 else if( i == 0 && j == 16 )
00232 {
00233
00234 strcpy( chLabel[k], "H2+ " );
00235 }
00236 else if( i == 0 && j == 19 )
00237 {
00238
00239 strcpy( chLabel[k], "Comp" );
00240 }
00241 else if( i == 0 && j == 22 )
00242 {
00243
00244 strcpy( chLabel[k], "line" );
00245 }
00246 else if( i == 0 && j == 23 )
00247 {
00248
00249
00250 strcpy( chLabel[k], "Hlin" );
00251 }
00252 else if( i == 0 && j == 24 )
00253 {
00254
00255 strcpy( chLabel[k], "ChaT" );
00256 }
00257 else if( i == 1 && j == 3 )
00258 {
00259
00260 strcpy( chLabel[k], "He3l" );
00261 }
00262 else if( i == 1 && j == 5 )
00263 {
00264
00265 strcpy( chLabel[k], "adve" );
00266 }
00267 else if( i == 1 && j == 6 )
00268 {
00269
00270 strcpy( chLabel[k], "CR H" );
00271 }
00272 else if( i == 25 && j == 27 )
00273 {
00274
00275 strcpy( chLabel[k], "Fe 2" );
00276 }
00277 else
00278 {
00279 sprintf( chLabel[k], "[%ld][%ld]" , i , j );
00280 }
00281 }
00282
00283
00284
00285 spsort(
00286
00287 save,
00288
00289 ipnt,
00290
00291 ipOrdered,
00292
00293
00294 -1,
00295
00296 &nFail);
00297
00298
00299
00300
00301 fprintf( io, "%.5e\t%.4e\t%.4e\t%.4e",
00302 radius.depth_mid_zone,
00303 phycon.te,
00304 heat_total,
00305 cool_total );
00306
00307
00308 ipnt = MIN2( ipnt , IPRINT );
00309
00310 for( k=0; k < ipnt; k++ )
00311 {
00312 int ip = ipOrdered[k];
00313 i = ipsave[ip];
00314 j = jpsave[ip];
00315 ASSERT( i<LIMELM && j<LIMELM );
00316
00317 if(k > 4 && thermal.heating[i][j]/heat_total < punch.WeakHeatCool)
00318 break;
00319
00320 fprintf( io, "\t%s\t%.7f ",
00321 chLabel[ip], save[ip] );
00322 }
00323 fprintf( io, " \n" );
00324
00325
00326
00327 lgHeatLine = false;
00328
00329
00330 for( i=0; i < ipnt; i++ )
00331 {
00332
00333 if( ipsave[i] == 0 && jpsave[i] == 22 )
00334 lgHeatLine = true;
00335 }
00336
00337 if( lgHeatLine )
00338 {
00339
00340 FndLineHt(&level,&ipStrong,&Strong);
00341 if( Strong/heat_total > 0.005 )
00342 {
00343 if( level == 1 )
00344 {
00345 strcpy( chLbl, chLineLbl(&TauLines[ipStrong] ) );
00346 TauIn = TauLines[ipStrong].TauIn;
00347 Pump = TauLines[ipStrong].pump;
00348 EscP = TauLines[ipStrong].Pesc;
00349 CS = TauLines[ipStrong].cs;
00350
00351 ColHeat = TauLines[ipStrong].heat/
00352 heat_total;
00353 }
00354 else if( level == 2 )
00355 {
00356 strcpy( chLbl, chLineLbl(&TauLine2[ipStrong]) );
00357 TauIn = TauLine2[ipStrong].TauIn;
00358 Pump = TauLine2[ipStrong].pump;
00359 EscP = TauLine2[ipStrong].Pesc;
00360 CS = TauLine2[ipStrong].cs;
00361 ColHeat = TauLine2[ipStrong].heat/
00362 heat_total;
00363 }
00364 else if( level == 3 )
00365
00366 {
00367 strcpy( chLbl, chLineLbl(&C12O16Rotate[ipStrong]) );
00368 TauIn = C12O16Rotate[ipStrong].TauIn;
00369 Pump = C12O16Rotate[ipStrong].pump;
00370 EscP = C12O16Rotate[ipStrong].Pesc;
00371 CS = C12O16Rotate[ipStrong].cs;
00372 ColHeat = C12O16Rotate[ipStrong].heat/
00373 heat_total;
00374 }
00375 else if( level == 4 )
00376
00377 {
00378 strcpy( chLbl, chLineLbl(&C13O16Rotate[ipStrong]) );
00379 TauIn = C13O16Rotate[ipStrong].TauIn;
00380 Pump = C13O16Rotate[ipStrong].pump;
00381 EscP = C13O16Rotate[ipStrong].Pesc;
00382 CS = C13O16Rotate[ipStrong].cs;
00383 ColHeat = C13O16Rotate[ipStrong].heat/
00384 heat_total;
00385 }
00386 else if( level == 5 )
00387
00388 {
00389 strcpy( chLbl, chLineLbl(&HFLines[ipStrong]) );
00390 TauIn = HFLines[ipStrong].TauIn;
00391 Pump = HFLines[ipStrong].pump;
00392 EscP = HFLines[ipStrong].Pesc;
00393 CS = HFLines[ipStrong].cs;
00394 ColHeat = HFLines[ipStrong].heat/
00395 heat_total;
00396 }
00397 else
00398 {
00399 fprintf( ioQQQ, " HeatPunch insane level%4ld\n",
00400 level );
00401 puts( "[Stop in HeatPunch]" );
00402 cdEXIT(EXIT_FAILURE);
00403 }
00404 fprintf( io, " LHeat lv%2ld %10.10s TIn%10.2e Pmp%9.1e EscP%9.1e CS%9.1e Hlin/tot%10.2e\n",
00405 level, chLbl, TauIn, Pump, EscP, CS, ColHeat );
00406 }
00407 }
00408 for( i=0; i<LIMELM*LIMELM; ++i )
00409 {
00410 free(chLabel[i]);
00411 }
00412
00413 free(chLabel);
00414 free(ipOrdered);
00415 free(jpsave);
00416 free(ipsave);
00417 free(save);
00418
00419 DEBUG_EXIT( "HeatPunch()" );
00420 return;
00421 }