00001
00002
00003
00004 #include "cddefines.h"
00005 #include "thermal.h"
00006 #include "dynamics.h"
00007 #include "radius.h"
00008 #include "conv.h"
00009 #include "phycon.h"
00010 #include "save.h"
00011 #include "grainvar.h"
00012 #include "hmi.h"
00013 #include "coolheavy.h"
00014
00015
00016
00017 static const int IPRINT = 100;
00018
00019
00020 void CoolSave(FILE * io, char chJob[])
00021 {
00022 long int i,
00023 ip,
00024 is;
00025
00026 int nFail;
00027
00028 double cset,
00029 cool_total,
00030 heat_total;
00031
00032 realnum
00033 *csav,
00034 *sgnsav;
00035 long int *index;
00036
00037 DEBUG_ENTRY( "CoolSave()" );
00038
00039
00040 index = (long int *)CALLOC((size_t)thermal.ncltot,sizeof(long int));
00041 csav = (realnum *)CALLOC((size_t)thermal.ncltot,sizeof(realnum));
00042 sgnsav = (realnum *)CALLOC((size_t)thermal.ncltot,sizeof(realnum));
00043
00044 cool_total = thermal.ctot;
00045 heat_total = thermal.htot;
00046
00047
00048
00049
00050
00051
00052 cool_total -= dynamics.Cool();
00053 heat_total -= dynamics.Heat();
00054 # if 0
00055 if(dynamics.Cool > dynamics.Heat())
00056 {
00057 cool_total -= dynamics.Heat();
00058 heat_total -= dynamics.Heat();
00059 }
00060 else
00061 {
00062 cool_total -= dynamics.Cool;
00063 heat_total -= dynamics.Cool;
00064 }
00065 # endif
00066
00067
00068
00069
00070 cset = cool_total*save.WeakHeatCool;
00071
00072
00073 ip = thermal.ncltot;
00074
00075 for( i=0; i < ip; i++ )
00076 {
00077 csav[i] = (realnum)(MAX2(thermal.cooling[i],thermal.heatnt[i])/
00078 SDIV(cool_total));
00079
00080
00081 if( thermal.heatnt[i] == 0. )
00082 {
00083 sgnsav[i] = 1.;
00084 }
00085 else
00086 {
00087 sgnsav[i] = -1.;
00088 }
00089 }
00090
00091
00092
00093
00094 spsort(
00095
00096 csav,
00097
00098 ip,
00099
00100 index,
00101
00102
00103 -1,
00104
00105 &nFail);
00106
00107
00108 if( !conv.lgConvTemp )
00109 {
00110 fprintf( io, "#>>>> Temperature not converged.\n" );
00111 }
00112 else if( !conv.lgConvEden )
00113 {
00114 fprintf( io, "#>>>> Electron density not converged.\n" );
00115 }
00116 else if( !conv.lgConvIoniz() )
00117 {
00118 fprintf( io, "#>>>> Ionization not converged.\n" );
00119 }
00120 else if( !conv.lgConvPres )
00121 {
00122 fprintf( io, "#>>>> Pressure not converged.\n" );
00123 }
00124
00125 if( strcmp(chJob,"EACH") == 0 )
00126 {
00127
00128 fprintf( io, "%.5e\t%.4e\t%.4e",
00129 radius.depth_mid_zone,
00130 phycon.te,
00131 cool_total );
00132 double debug_ctot = 0.;
00133
00134 for( int i=0 ; i <= LIMELM ; i++)
00135 {
00136 fprintf( io, "\t%.4e", thermal.elementcool[i] );
00137 debug_ctot += thermal.elementcool[i];
00138 }
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158 fprintf( io, "\t%.4e", MAX2(0.,gv.GasCoolColl) );
00159 debug_ctot += MAX2(0.,gv.GasCoolColl);
00160
00161 fprintf( io, "\t%.4e", MAX2(0.,-hmi.HeatH2Dexc_used) );
00162 debug_ctot += MAX2(0.,-hmi.HeatH2Dexc_used);
00163
00164 fprintf( io, "\t%.4e", thermal.char_tran_cool );
00165 debug_ctot += thermal.char_tran_cool;
00166
00167 fprintf( io, "\t%.4e", hmi.hmicol );
00168 debug_ctot += hmi.hmicol;
00169
00170 fprintf( io, "\t%.4e", CoolHeavy.h2line );
00171 debug_ctot += CoolHeavy.h2line;
00172
00173 fprintf( io, "\t%.4e", CoolHeavy.HD );
00174 debug_ctot += CoolHeavy.HD;
00175
00176 fprintf( io, "\t%.4e", CoolHeavy.H2PlsCool );
00177 debug_ctot += CoolHeavy.H2PlsCool;
00178
00179 fprintf( io, "\t%.4e", MAX2(0.,CoolHeavy.brems_cool_net) );
00180 debug_ctot += MAX2(0.,CoolHeavy.brems_cool_net);
00181
00182 fprintf( io, "\t%.4e", CoolHeavy.heavfb );
00183 debug_ctot += CoolHeavy.heavfb;
00184
00185 fprintf( io, "\t%.4e", CoolHeavy.eebrm );
00186 debug_ctot += CoolHeavy.eebrm;
00187
00188 fprintf( io, "\t%.4e", CoolHeavy.tccool );
00189 debug_ctot += CoolHeavy.tccool;
00190
00191 fprintf( io, "\t%.4e", CoolHeavy.cextxx );
00192 debug_ctot += CoolHeavy.cextxx;
00193
00194 fprintf( io, "\t%.4e", CoolHeavy.expans );
00195 debug_ctot += CoolHeavy.expans;
00196
00197 fprintf( io, "\t%.4e", CoolHeavy.cyntrn );
00198 debug_ctot += CoolHeavy.cyntrn;
00199
00200 fprintf( io, "\t%.4e", CoolHeavy.colmet );
00201 debug_ctot += CoolHeavy.colmet;
00202
00203 fprintf( io, "\t%.4e", thermal.dima );
00204 debug_ctot += thermal.dima;
00205 fprintf( io, " \n" );
00206 if( fabs( (debug_ctot - cool_total)/cool_total ) > 1e-10 )
00207 {
00208 fprintf( ioQQQ , "PROBLEM with the SAVE EACH COOLING output\n" );
00209 fprintf( ioQQQ , "PROBLEM One or more coolants have been lost, the sum of the reported cooling is %.4e\n", debug_ctot );
00210 fprintf( ioQQQ , "PROBLEM The total cooling is %.4ee\n", cool_total );
00211 fprintf( ioQQQ , "PROBLEM The difference is %.4e\n", cool_total - debug_ctot );
00212 cdEXIT(EXIT_FAILURE);
00213 }
00214 }
00215 else if( strcmp(chJob,"COOL") == 0 )
00216 {
00217
00218
00219
00220 fprintf( io, "%.5e\t%.4e\t%.4e\t%.4e",
00221 radius.depth_mid_zone,
00222 phycon.te,
00223 heat_total,
00224 cool_total );
00225
00226
00227 ip = MIN2( ip , IPRINT );
00228
00229
00230
00231
00232 for( is=0; is < ip; is++ )
00233 {
00234 if(is > 4 && (thermal.cooling[index[is]] < cset && thermal.heatnt[index[is]] < cset))
00235 break;
00236 fprintf( io, "\t%s %.1f\t%.7f",
00237 thermal.chClntLab[index[is]],
00238 thermal.collam[index[is]],
00239 sign(csav[index[is]],sgnsav[index[is]]) );
00240 }
00241 fprintf( io, " \n" );
00242 }
00243 else
00244 TotalInsanity();
00245
00246
00247 free(sgnsav);
00248 free(csav);
00249 free(index);
00250
00251 return;
00252 }
00253