00001
00002
00003
00004 #include "cddefines.h"
00005 #include "physconst.h"
00006 #include "taulines.h"
00007 #include "iso.h"
00008 #include "secondaries.h"
00009 #include "conv.h"
00010 #include "elementnames.h"
00011 #include "atmdat.h"
00012 #include "rfield.h"
00013 #include "dense.h"
00014 #include "trace.h"
00015 #include "thirdparty.h"
00016 #include "hydrogenic.h"
00017 #include "dynamics.h"
00018
00019
00020 static double **z;
00021
00022 static void PrtHydroTrace2( long int nelem )
00023 {
00024 long ipHi , j , n;
00025 double collider;
00026
00027 DEBUG_ENTRY( "PrtHydroTrace2()" );
00028
00029
00030
00031
00032 if( nelem==ipHYDROGEN )
00033 {
00034
00035 collider = dense.EdenHontoHCorr;
00036 }
00037 else
00038 {
00039 collider = dense.EdenHCorr;
00040 }
00041
00042 if( nelem == ipHYDROGEN )
00043 {
00044 double sum=0 , HRateDestGnd[10];
00045
00046 for(ipHi=1; ipHi<iso.numLevels_local[ipH_LIKE][nelem]; ++ipHi )
00047 {
00048 sum += EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].pump;
00049 }
00050 fprintf(ioQQQ,"DEBUG pump\t%.2f", fnzone);
00051 for(ipHi=1; ipHi<iso.numLevels_local[ipH_LIKE][nelem]; ++ipHi )
00052 {
00053 fprintf(ioQQQ,"\t%.3f",EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].pump/SDIV(sum));
00054 }
00055 fprintf(ioQQQ,"\n");
00056 ipHi = iso.numLevels_local[ipH_LIKE][nelem]-2;
00057 fprintf(ioQQQ,"DEBUG pumphi\t%.2f\t%.3e\t%.3e\t%.3e\n",
00058 fnzone,
00059 EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].pump,
00060 rfield.OccNumbIncidCont[EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].ipCont-1],
00061 rfield.OccNumbDiffCont[EmisLines[ipH_LIKE][nelem][ipHi][ipH1s].ipCont-1]);
00062
00063 HRateDestGnd[0] = iso.RateLevel2Cont[ipH_LIKE][nelem][ipH1s];
00064 HRateDestGnd[1] = iso.gamnc[ipH_LIKE][nelem][ipH1s]/
00065 iso.RateLevel2Cont[ipH_LIKE][nelem][ipH1s];
00066
00067 HRateDestGnd[2] = iso.ColIoniz[ipH_LIKE][nelem][ipH1s]*
00068 collider/iso.RateLevel2Cont[ipH_LIKE][nelem][ipH1s];
00069 HRateDestGnd[3] = secondaries.csupra[nelem][nelem]/iso.RateLevel2Cont[ipH_LIKE][nelem][ipH1s];
00070 HRateDestGnd[4] = secondaries.Hx12[MIN2(nelem,1)][ipH2p]*
00071 9./iso.RateLevel2Cont[ipH_LIKE][nelem][ipH1s];
00072 HRateDestGnd[5] = atmdat.HCharExcIonTotal/
00073 iso.RateLevel2Cont[ipH_LIKE][nelem][ipH1s];
00074 HRateDestGnd[6] = atmdat.HCharExcIonTotal;
00075 HRateDestGnd[7] = sum/iso.RateLevel2Cont[ipH_LIKE][nelem][ipH1s];
00076 HRateDestGnd[8] = sum;
00077
00078 fprintf(ioQQQ," grnd dest fracs\t%.2f",fnzone);
00079 for( j=ipH1s; j < 9; j++ )
00080 {
00081 fprintf( ioQQQ,"\t");
00082 fprintf( ioQQQ,PrintEfmt("%8.1e", HRateDestGnd[j] ));
00083 }
00084 fprintf(ioQQQ,"\thcoldc\t%e\n", dense.EdenHCorr );
00085 fprintf( ioQQQ, "\n" );
00086 }
00087
00088 fprintf( ioQQQ, " pop level others => (HydroLevelPop)\n" );
00089
00090 for( n=ipH1s; n < iso.numLevels_local[ipH_LIKE][nelem]; n++ )
00091 {
00092 fprintf( ioQQQ, " HII%2ld", n );
00093 for( j=ipH1s; j < iso.numLevels_local[ipH_LIKE][nelem]; j++ )
00094 {
00095 fprintf( ioQQQ," ");
00096
00097 fprintf( ioQQQ,"%.9e\t", z[j][n] );
00098 }
00099 fprintf( ioQQQ, "\n" );
00100 }
00101
00102 DEBUG_EXIT( "HydroLevelPop()" );
00103 return;
00104 }
00105
00106 #if defined (__ICC) && defined(__ia64) && __INTEL_COMPILER < 910
00107 #pragma optimize("", off)
00108 #endif
00109
00110 void HydroLevelPop(long int nelem)
00111 {
00112 long int n,
00113 ipHi,
00114 ipLo,
00115 i,
00116 j,
00117 level,
00118 level_error;
00119 double *amat;
00120 double collider;
00121
00122 double BigError;
00123
00124 int32 nerror;
00125 int32 *ipiv;
00126
00127 double
00128 *creation,
00129 *Save_creation,
00130 **SaveZ,
00131 *work;
00132 double *error;
00133
00134 DEBUG_ENTRY( "HydroLevelPop()" );
00135
00136
00137 ASSERT( nelem >= 0 );
00138 ASSERT( nelem < LIMELM );
00139
00140 ipiv = (int32 *)MALLOC(sizeof(int32)*(unsigned)(iso.numLevels_max[ipH_LIKE][nelem]) );
00141 creation = (double *)MALLOC(sizeof(double)*(unsigned)(iso.numLevels_max[ipH_LIKE][nelem]) );
00142 Save_creation = (double *)MALLOC(sizeof(double)*(unsigned)(iso.numLevels_max[ipH_LIKE][nelem]) );
00143 work = (double *)MALLOC(sizeof(double)*(unsigned)(iso.numLevels_max[ipH_LIKE][nelem]) );
00144 error = (double *)MALLOC(sizeof(double)*(unsigned)(iso.numLevels_max[ipH_LIKE][nelem]) );
00145
00146
00147 SaveZ = (double **)MALLOC(sizeof(double *)*(unsigned)(iso.numLevels_max[ipH_LIKE][nelem]) );
00148
00149 z = (double **)MALLOC(sizeof(double *)*(unsigned)(iso.numLevels_max[ipH_LIKE][nelem]) );
00150
00151
00152 for( i=0; i<(iso.numLevels_max[ipH_LIKE][nelem]); ++i )
00153 {
00154 SaveZ[i] = (double *)MALLOC(sizeof(double)*(unsigned)(iso.numLevels_max[ipH_LIKE][nelem]) );
00155
00156 z[i] = (double *)MALLOC(sizeof(double)*(unsigned)(iso.numLevels_max[ipH_LIKE][nelem]) );
00157 }
00158
00159
00160
00161
00162 if( nelem==ipHYDROGEN )
00163 {
00164
00165 collider = dense.EdenHontoHCorr;
00166 }
00167 else
00168 {
00169 collider = dense.EdenHCorr;
00170 }
00171
00172
00173
00174 for( level=ipH1s; level < iso.numLevels_local[ipH_LIKE][nelem]; ++level )
00175 {
00176
00177 creation[level] = iso.RateCont2Level[ipH_LIKE][nelem][level];
00178 }
00179
00180
00181
00182 for( level=ipH1s; level < iso.numLevels_local[ipH_LIKE][nelem]; level++ )
00183 {
00184
00185
00186 z[level][level] = iso.RateLevel2Cont[ipH_LIKE][nelem][level];
00187
00188
00189 for( ipLo=ipH1s; ipLo < level; ipLo++ )
00190 {
00191
00192
00193 double CollRate = EmisLines[ipH_LIKE][nelem][level][ipLo].ColUL* collider;
00194
00195 z[ipLo][level] =
00196 -CollRate *
00197 (double)iso.stat[ipH_LIKE][nelem][level]/(double)iso.stat[ipH_LIKE][nelem][ipLo]*
00198 iso.Boltzmann[ipH_LIKE][nelem][level][ipLo] -
00199 EmisLines[ipH_LIKE][nelem][level][ipLo].pump *
00200 KillIfBelowPlasma(EmisLines[ipH_LIKE][nelem][level][ipLo].EnergyWN/WAVNRYD);
00201
00202
00203 z[level][level] += (double)EmisLines[ipH_LIKE][nelem][level][ipLo].pump *
00204 (double)iso.stat[ipH_LIKE][nelem][ipLo]/(double)iso.stat[ipH_LIKE][nelem][level] *
00205 KillIfBelowPlasma(EmisLines[ipH_LIKE][nelem][level][ipLo].EnergyWN/WAVNRYD);
00206
00207
00208 z[level][level] += CollRate;
00209
00210
00211 z[level][level] +=
00212 EmisLines[ipH_LIKE][nelem][level][ipLo].Aul*
00213 (EmisLines[ipH_LIKE][nelem][level][ipLo].Pesc +
00214 EmisLines[ipH_LIKE][nelem][level][ipLo].Pelec_esc +
00215 EmisLines[ipH_LIKE][nelem][level][ipLo].Pdest) *
00216 KillIfBelowPlasma(EmisLines[ipH_LIKE][nelem][level][ipLo].EnergyWN/WAVNRYD);
00217 }
00218
00219
00220
00221 for( ipHi=level + 1; ipHi < iso.numLevels_local[ipH_LIKE][nelem]; ipHi++ )
00222 {
00223 double RadDecay , CollRate;
00224 RadDecay =
00225 EmisLines[ipH_LIKE][nelem][ipHi][level].Aul*
00226 (EmisLines[ipH_LIKE][nelem][ipHi][level].Pesc +
00227 EmisLines[ipH_LIKE][nelem][ipHi][level].Pelec_esc +
00228 EmisLines[ipH_LIKE][nelem][ipHi][level].Pdest);
00229
00230
00231
00232
00233
00234 CollRate = EmisLines[ipH_LIKE][nelem][ipHi][level].ColUL *collider;
00235
00236 z[ipHi][level] =
00237 -( RadDecay +
00238 (double)EmisLines[ipH_LIKE][nelem][ipHi][level].pump *
00239 (double)iso.stat[ipH_LIKE][nelem][level]/(double)iso.stat[ipH_LIKE][nelem][ipHi] ) *
00240 KillIfBelowPlasma(EmisLines[ipH_LIKE][nelem][ipHi][level].EnergyWN/WAVNRYD) -
00241 CollRate;
00242
00243
00244 z[level][level] += EmisLines[ipH_LIKE][nelem][ipHi][level].pump *
00245 KillIfBelowPlasma(EmisLines[ipH_LIKE][nelem][ipHi][level].EnergyWN/WAVNRYD);
00246
00247
00248 z[level][level] += (double)iso.stat[ipH_LIKE][nelem][ipHi] / (double)iso.stat[ipH_LIKE][nelem][level] *
00249 iso.Boltzmann[ipH_LIKE][nelem][ipHi][level]*CollRate;
00250 }
00251 }
00252
00253
00254
00255
00256
00257 z[ipH2s][ipH1s] -= iso.TwoNu_induc_dn[ipH_LIKE][nelem]*iso.lgInd2nu_On;
00258 z[ipH1s][ipH2s] -= iso.TwoNu_induc_up[ipH_LIKE][nelem]*iso.lgInd2nu_On;
00259
00260
00261 z[ipH1s][ipH1s] += iso.TwoNu_induc_up[ipH_LIKE][nelem]*iso.lgInd2nu_On;
00262 z[ipH2s][ipH2s] += iso.TwoNu_induc_dn[ipH_LIKE][nelem]*iso.lgInd2nu_On;
00263
00264
00265 if( secondaries.Hx12[0][1] * iso.lgColl_excite[ipH_LIKE] > 0. )
00266 {
00267
00268
00269 for( level=ipH2s; level < iso.numLevels_local[ipH_LIKE][nelem]; level++ )
00270 {
00271 double RateUp , RateDown;
00272
00273 RateUp = secondaries.Hx12[MIN2(nelem,1)][level];
00274
00275 RateDown = RateUp * (double)iso.stat[ipH_LIKE][nelem][ipH1s] /
00276 (double)iso.stat[ipH_LIKE][nelem][level];
00277
00278
00279
00280 z[ipH1s][ipH1s] += RateUp;
00281
00282
00283 z[level][ipH1s] -= RateDown;
00284
00285
00286 z[ipH1s][level] -= RateUp;
00287
00288 z[level][level] += RateDown;
00289 }
00290 }
00291
00292
00293
00294 {
00295
00296 enum {DEBUG_LOC=false};
00297
00298 if( DEBUG_LOC && nelem==ipHYDROGEN && nzone > 600 )
00299 {
00300 fprintf(ioQQQ, "DEBUG Hn=2\t%.2f\t%.3e\t%.3e\t%.3e\n",
00301 fnzone,
00302 z[ipH1s][ipH1s],
00303 z[ipH2p][ipH2p],
00304 EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].Pdest );
00305 }
00306 }
00307
00308 #ifdef FOO
00309 {
00310 float mytot,err;
00311 int myeq;
00312
00313
00314 for(myeq=ipH1s; myeq<iso.numLevels_local[ipH_LIKE][nelem]; myeq++ ) {
00315 mytot = 0.;
00316 for(level=ipH1s; level<iso.numLevels_local[ipH_LIKE][nelem]; level++ ) {
00317 mytot += z[myeq][level];
00318 }
00319 err = fabs((mytot-iso.RateLevel2Cont[ipH_LIKE][nelem][myeq])/z[myeq][myeq]);
00320 if(err > 1e-6)
00321 printf("BalChk: %d %d e %g\n",nelem,myeq,err);
00322 }
00323 }
00324 #endif
00325
00326
00327
00328 if( dynamics.lgAdvection )
00329 {
00330
00331
00332 creation[ipH1s] +=
00333 dynamics.Source[nelem][nelem]/SDIV(dense.xIonDense[nelem][nelem+1])*
00334 dynamics.lgISO[ipH_LIKE];
00335
00336
00337
00338 for( ipLo=ipH1s; ipLo<iso.numLevels_local[ipH_LIKE][nelem]; ++ipLo )
00339 {
00340 z[ipLo][ipLo] += dynamics.Rate*dynamics.lgISO[ipH_LIKE];
00341 }
00342 }
00343
00344
00345
00346
00347
00348
00349
00350 if( (trace.lgTrace && trace.lgIsoTraceFull[ipH_LIKE]) && (nelem == trace.ipIsoTrace[ipH_LIKE]) )
00351 PrtHydroTrace2(nelem);
00352
00353
00354
00355
00356 for( j=ipH1s; j < iso.numLevels_local[ipH_LIKE][nelem]; j++ )
00357 {
00358 for( n=ipH1s; n < iso.numLevels_local[ipH_LIKE][nelem]; n++ )
00359 {
00360 SaveZ[n][j] = z[n][j];
00361 }
00362 Save_creation[j] = creation[j];
00363 }
00364
00365
00366
00367 # ifdef AMAT
00368 # undef AMAT
00369 # endif
00370
00371 # define AMAT(I_,J_) (*(amat+(I_)*iso.numLevels_local[ipH_LIKE][nelem]+(J_)))
00372
00373
00374 amat=(double*)MALLOC( (sizeof(double)*(unsigned)(iso.numLevels_local[ipH_LIKE][nelem]*(iso.numLevels_local[ipH_LIKE][nelem])) ));
00375
00376
00377
00378 for( j=ipH1s; j < iso.numLevels_local[ipH_LIKE][nelem]; j++ )
00379 {
00380 for( n=ipH1s; n < iso.numLevels_local[ipH_LIKE][nelem]; n++ )
00381 {
00382
00383 AMAT(n,j) = z[n][j];
00384 }
00385 }
00386
00387
00388
00389
00390
00391
00392
00393
00394 nerror = 0;
00395
00396
00397 getrf_wrapper(iso.numLevels_local[ipH_LIKE][nelem],iso.numLevels_local[ipH_LIKE][nelem],
00398 amat,iso.numLevels_local[ipH_LIKE][nelem],ipiv,&nerror);
00399
00400 getrs_wrapper('N',iso.numLevels_local[ipH_LIKE][nelem],1,amat,iso.numLevels_local[ipH_LIKE][nelem],ipiv,
00401 creation,iso.numLevels_local[ipH_LIKE][nelem],&nerror);
00402
00403
00404
00405
00406
00407
00408
00409
00410 if( nerror != 0 )
00411 {
00412 fprintf( ioQQQ, " HydroLevelPop: dgetrs finds singular or ill-conditioned matrix\n" );
00413 puts( "[Stop in HydroLevelPop]" );
00414 cdEXIT(EXIT_FAILURE);
00415 }
00416
00417 free(amat);
00418
00419
00420
00421 for( level=ipH1s; level < iso.numLevels_local[ipH_LIKE][nelem]; level++ )
00422 {
00423 double BigSoln= 0.;
00424 error[level] = 0.;
00425 for( n=ipH1s; n < iso.numLevels_local[ipH_LIKE][nelem]; n++ )
00426 {
00427
00428 if( fabs(SaveZ[n][level] ) > BigSoln )
00429 BigSoln = fabs(SaveZ[n][level]);
00430
00431 error[level] += SaveZ[n][level]*creation[n];
00432 }
00433
00434 if( BigSoln > 0. )
00435 {
00436 error[level] = (error[level] - Save_creation[level])/ BigSoln;
00437 }
00438 else
00439 {
00440 error[level] = 0.;
00441 }
00442 }
00443
00444
00445 BigError = -1.;
00446 level_error = -1;
00447
00448 for( level=ipH1s; level < iso.numLevels_local[ipH_LIKE][nelem]; level++ )
00449 {
00450 double abserror;
00451 abserror = fabs( error[level]);
00452
00453 if( abserror > BigError )
00454 {
00455 BigError = abserror;
00456 level_error = level;
00457 }
00458 }
00459
00460
00461
00462 if( BigError > FLT_EPSILON )
00463 {
00464 if( !conv.lgSearch )
00465 fprintf(ioQQQ," PROBLEM " );
00466
00467 fprintf(ioQQQ,
00468 " HydroLevelPop zone %.2f - largest residual in hydrogenic %s nelem=%li matrix inversion is %g "
00469 "level was %li Search?%c \n",
00470 fnzone,
00471 elementnames.chElementName[nelem],
00472 nelem ,
00473 BigError ,
00474 level_error,
00475 TorF(conv.lgSearch) );
00476
00477
00478
00479
00480
00481 }
00482
00483
00484
00485
00486
00487 for( ipLo=ipH1s; ipLo < iso.numLevels_local[ipH_LIKE][nelem]; ipLo++ )
00488 {
00489 iso.Pop2Ion[ipH_LIKE][nelem][ipLo] = creation[ipLo];
00490 if( iso.Pop2Ion[ipH_LIKE][nelem][ipLo] <= 0 )
00491 {
00492 fprintf(ioQQQ," non-positive level pop for H iso, nelem = %li = %s, level=%li val=%.3e\n",
00493 nelem ,
00494 elementnames.chElementSym[nelem],
00495 ipLo,
00496 iso.Pop2Ion[ipH_LIKE][nelem][ipLo] );
00497 }
00498
00499 if( iso.PopLTE[ipH_LIKE][nelem][ipLo] > 0. )
00500 {
00501 iso.DepartCoef[ipH_LIKE][nelem][ipLo] =
00502 (iso.Pop2Ion[ipH_LIKE][nelem][ipLo]/
00503 (iso.PopLTE[ipH_LIKE][nelem][ipLo]* dense.eden) );
00504 }
00505 else
00506 {
00507 iso.DepartCoef[ipH_LIKE][nelem][ipLo] = 0.;
00508 }
00509 }
00510
00511 free( ipiv );
00512 free( work );
00513 free( creation );
00514 free( Save_creation );
00515 free( error );
00516
00517
00518
00519 for( i=0; i<iso.numLevels_max[ipH_LIKE][nelem]; ++i )
00520 {
00521 free( SaveZ[i] );
00522 free( z[i] );
00523 }
00524 free( SaveZ );
00525 free( z );
00526
00527 DEBUG_EXIT( "HydroLevelPop()" );
00528 return;
00529 }
00530
00531