00001
00002
00003
00004
00005 #include "cddefines.h"
00006 #include "phycon.h"
00007 #include "dense.h"
00008 #include "taulines.h"
00009 #include "input.h"
00010 #include "h2.h"
00011 #include "h2_priv.h"
00012 #include "mole.h"
00013
00014
00015 #define PRT_COLL false
00016
00017
00018
00019 #define N_H2_HE_FIT_PAR 8
00020 static realnum ***H2_He_coll_fit_par;
00021 static bool **lgDefn_H2He_coll;
00022
00023
00024
00025 #define N_H2_ORH2_FIT_PAR 8
00026 static realnum ***H2_ORH2_coll_fit_par;
00027 static bool **lgDefn_H2ORH2_coll;
00028
00029
00030
00031
00032 #define N_H2_PAH2_FIT_PAR 8
00033 static realnum ***H2_PAH2_coll_fit_par;
00034 static bool **lgDefn_H2PAH2_coll;
00035
00036
00037 STATIC realnum H2_CollidRateEvalOne(
00038
00039
00040 long iVibHi, long iRotHi,long iVibLo,
00041
00042 long iRotLo, long ipHi , long ipLo ,
00043
00044 long nColl )
00045 {
00046 double fitted;
00047 realnum rate;
00048 double t3Plus1 = phycon.te/1000. + 1.;
00049 double t3Plus1Squared = POW2(t3Plus1);
00050
00051
00052 double gbarcoll[N_X_COLLIDER][3] =
00053 {
00054 {-9.9265 , -0.1048 , 0.456 },
00055 {-8.281 , -0.1303 , 0.4931 },
00056 {-10.0357, -0.0243 , 0.67 },
00057 {-8.6213 , -0.1004 , 0.5291 },
00058 {-9.2719 , -0.0001 , 1.0391 }
00059 };
00060
00061 DEBUG_ENTRY( "H2_CollidRateEvalOne()" );
00062
00063
00064
00065 if( nColl == 1 && mole.lgH2_He_ORNL )
00066 {
00067
00068
00069
00070
00071
00072
00073
00074
00075 if( (fitted=H2_He_coll(ipHi , ipLo , phycon.te ))>0. )
00076 {
00077
00078
00079
00080
00081 rate = (realnum)fitted*mole.lgColl_deexec_Calc;
00082
00083
00084
00085
00086
00087
00088 }
00089
00090
00091
00092
00093 else if( mole.lgColl_gbar &&
00094 (H2_lgOrtho[0][iVibHi][iRotHi]-H2_lgOrtho[0][iVibLo][iRotLo]==0) )
00095 {
00096
00097
00098 double ediff = energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo];
00099
00100
00101
00102 ediff = MAX2(100., ediff );
00103 rate = (realnum)pow(10. ,
00104 gbarcoll[nColl][0] + gbarcoll[nColl][1] *
00105 pow(ediff,gbarcoll[nColl][2]) )*mole.lgColl_deexec_Calc;
00106
00107
00108
00109
00110
00111
00112
00113 }
00114 else
00115 rate = 0;
00116 }
00117 else if( nColl==4 )
00118 {
00119
00120
00121
00122 if( CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][1] != 0 )
00123 {
00124 rate = CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][0] * 1e-10f *
00125
00126 (realnum)sexp( CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][1]/1000./phycon.te_eV)*mole.lgColl_deexec_Calc;
00127
00128
00129
00130
00131
00132
00133 }
00134
00135
00136
00137 else if( mole.lgColl_gbar )
00138 {
00139
00140
00141 double ediff = energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo];
00142 ediff = MAX2(100., ediff );
00143 rate = (realnum)pow(10. ,
00144 gbarcoll[nColl][0] + gbarcoll[nColl][1] *
00145 pow(ediff ,gbarcoll[nColl][2]) )*mole.lgColl_deexec_Calc;
00146
00147
00148
00149
00150
00151
00152
00153 }
00154 else
00155 rate = 0;
00156 }
00157
00158 else if( nColl == 2 && mole.lgH2_ORH2_ORNL )
00159 {
00160
00161
00162
00163
00164
00165 if( (fitted=H2_ORH2_coll(ipHi , ipLo , phycon.te ))>0. )
00166 {
00167 rate = (realnum)fitted*mole.lgColl_deexec_Calc;
00168 if( PRT_COLL )
00169 fprintf(ioQQQ,"col fit\t%li\t%li\t%.4e\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
00170 ipLo,ipHi,phycon.te,nColl,
00171 iVibHi,iRotHi,iVibLo,iRotLo,
00172 energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
00173 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] );
00174 }
00175 else
00176 rate = 0;
00177 }
00178
00179
00180 else if( nColl == 3 && mole.lgH2_PAH2_ORNL )
00181 {
00182
00183
00184
00185
00186
00187 if( (fitted=H2_PAH2_coll(ipHi , ipLo , phycon.te ))>0. )
00188 {
00189 rate = (realnum)fitted*mole.lgColl_deexec_Calc;
00190 if( PRT_COLL )
00191 fprintf(ioQQQ,"col fit\t%li\t%li\t%.4e\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
00192 ipLo,ipHi,phycon.te,nColl,
00193 iVibHi,iRotHi,iVibLo,iRotLo,
00194 energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
00195 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] );
00196 }
00197 else
00198 rate = 0;
00199 }
00200
00201 else if( CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][0]!= 0 )
00202 {
00203
00204
00205
00206
00207
00208 double r = CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][0] +
00209 CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][1]/t3Plus1 +
00210 CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][2]/t3Plus1Squared;
00211
00212 rate = (realnum)pow(10.,r)*mole.lgColl_deexec_Calc;
00213
00214
00215
00216
00217
00218
00219
00220 }
00221
00222
00223
00224 else if( mole.lgColl_gbar &&
00225 (H2_lgOrtho[0][iVibHi][iRotHi]-H2_lgOrtho[0][iVibLo][iRotLo]==0) )
00226 {
00227
00228
00229 double ediff = energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo];
00230
00231
00232 ediff = MAX2(100., ediff );
00233 rate = (realnum)pow(10. ,
00234 gbarcoll[nColl][0] + gbarcoll[nColl][1] *
00235 pow(ediff,gbarcoll[nColl][2]) )*mole.lgColl_deexec_Calc;
00236
00237 if( nColl == 0 && h2.lgH2_H_coll_07 )
00238 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] *= 100.;
00239
00240
00241
00242
00243
00244
00245
00246 }
00247 else
00248 rate = 0;
00249
00250
00251
00252 if( !mole.lgH2_ortho_para_coll_on &&
00253 (H2_lgOrtho[0][iVibHi][iRotHi]-H2_lgOrtho[0][iVibLo][iRotLo]) )
00254 rate = 0.;
00255
00256 {
00257 enum {DEBUG_LOC=false};
00258 if( DEBUG_LOC )
00259 {
00260 fprintf(ioQQQ,"bugcoll\tiVibHi\t%li\tiRotHi\t%li\tiVibLo\t%li\tiRotLo\t%li\tcoll\t%.2e\n",
00261 iVibHi,iRotHi,iVibLo,iRotLo,
00262 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] );
00263 }
00264 }
00265 return rate;
00266 }
00267
00268
00269 void H2_CollidRateEvalAll( void )
00270 {
00271 long int numb_coll_trans = 0;
00272 double excit;
00273 long int iElecHi , iElecLo , ipHi , iVibHi , iRotHi ,
00274 ipLo , iVibLo , iRotLo , nColl;
00275
00276 DEBUG_ENTRY( "H2_CollidRateEvalAll()" );
00277
00278 if( PRT_COLL )
00279 fprintf(ioQQQ,"H2 coll deex rate coef\n"
00280 "VibHi\tRotHi\tVibLo\tRotLo\trate\n");
00281
00282 iElecHi = 0;
00283 iElecLo = 0;
00284 if(mole.nH2_TRACE >= mole.nH2_trace_full)
00285 fprintf(ioQQQ,"H2 set collision rates\n");
00286
00287
00288
00289
00290
00291 H2_coll_dissoc_rate_coef[0][0] = 0.;
00292 H2_coll_dissoc_rate_coef_H2[0][0] = 0.;
00293 for( ipHi=0; ipHi<nLevels_per_elec[0]; ++ipHi )
00294 {
00295 double energy;
00296
00297
00298 long int ip = H2_ipX_ener_sort[ipHi];
00299 iVibHi = ipVib_H2_energy_sort[ip];
00300 iRotHi = ipRot_H2_energy_sort[ip];
00301
00302
00303
00304
00305 energy = H2_DissocEnergies[0] - energy_wn[0][iVibHi][iRotHi];
00306 ASSERT( energy > 0. );
00307
00308 H2_coll_dissoc_rate_coef[iVibHi][iRotHi] =
00309 1e-14f * (realnum)sexp(energy/phycon.te_wn) * mole.lgColl_dissoc_coll;
00310
00311
00312
00313 H2_coll_dissoc_rate_coef_H2[iVibHi][iRotHi] =
00314 1e-11f * (realnum)sexp(energy/phycon.te_wn) * mole.lgColl_dissoc_coll;
00315
00316
00317
00318
00319
00320
00321 for( ipLo=0; ipLo<ipHi; ++ipLo )
00322 {
00323 ip = H2_ipX_ener_sort[ipLo];
00324 iVibLo = ipVib_H2_energy_sort[ip];
00325 iRotLo = ipRot_H2_energy_sort[ip];
00326
00327 ASSERT( energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo] > 0.);
00328
00329
00330
00331
00332
00333
00334
00335 ++numb_coll_trans;
00336
00337
00338
00339 if( PRT_COLL && iVibHi == 1 && iRotHi==3)
00340 fprintf(ioQQQ,"%li\t%li\t%li\t%li",
00341 iVibHi,iRotHi,iVibLo,iRotLo);
00342 for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
00343 {
00344 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] =
00345 H2_CollidRateEvalOne( iVibHi,iRotHi,iVibLo,iRotLo,
00346 ipHi , ipLo , nColl );
00347 if( PRT_COLL && iVibHi == 1 && iRotHi==3)
00348 fprintf(ioQQQ,"\t%.2e",H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] );
00349 }
00350 if( PRT_COLL && iVibHi == 1 && iRotHi==3)
00351 fprintf(ioQQQ,"\n");
00352 }
00353 }
00354 if( PRT_COLL )
00355 cdEXIT( EXIT_FAILURE );
00356
00357
00358
00359
00360
00361
00362 nColl = 0;
00363 iElecHi = 0;
00364 iElecLo = 0;
00365 iVibHi = 0;
00366 iVibLo = 0;
00367
00368
00369
00370 double excit1;
00371 double te_used = MAX2( 100. , phycon.te );
00372
00373 iRotLo = 0;
00374 iRotHi = 1;
00375 excit1 = sexp( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyK/te_used);
00376 excit = sexp( -(POW2(5.30-460./te_used)-21.2) )*1e-13;
00377
00378 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][0] = (realnum)(
00379 excit*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->g/
00380 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g /
00381
00382 SDIV(excit1) )*mole.lgColl_deexec_Calc *
00383
00384 mole.lgH2_ortho_para_coll_on;
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394 iRotLo = 0;
00395 iRotHi = 3;
00396 excit = sexp( -(POW2(6.36-373./te_used)-34.5) )*1e-13;
00397 excit1 = sexp( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyK/te_used);
00398 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][0] = (realnum)(
00399 excit*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->g/
00400 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g /
00401 SDIV(excit1) )*mole.lgColl_deexec_Calc *
00402
00403 mole.lgH2_ortho_para_coll_on;
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413 iRotLo = 1;
00414 iRotHi = 2;
00415 excit = sexp( -(POW2(5.35-454./te_used)-23.1 ) )*1e-13;
00416 excit1 = sexp( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyK/te_used);
00417 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][0] = (realnum)(
00418 excit*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->g/
00419 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g /
00420 SDIV(excit1) )*mole.lgColl_deexec_Calc *
00421
00422 mole.lgH2_ortho_para_coll_on;
00423
00424
00425
00426
00427
00428
00429
00430 if( phycon.te < 100. )
00431 {
00432
00433
00434 H2_CollRate[0][1][0][0][0] = (realnum)(H2_CollRate[0][0][1][0][0]*exp(-(3900-170.5)*(1./phycon.te - 1./100.)));
00435 H2_CollRate[0][3][0][0][0] = (realnum)(H2_CollRate[0][0][3][0][0]*exp(-(3900-1015.1)*(1./phycon.te - 1./100.)));
00436 H2_CollRate[0][2][0][1][0] = (realnum)(H2_CollRate[0][0][2][0][1]*exp(-(3900-339.3)*(1./phycon.te - 1./100.)));
00437 }
00438
00439
00440
00441
00442
00443
00444
00445
00446 if( mole.nH2_TRACE >= mole.nH2_trace_full )
00447 fprintf(ioQQQ,
00448 " collision rates updated for new temp, number of trans is %li\n",
00449 numb_coll_trans);
00450
00451
00452
00453
00454 return;
00455 }
00456
00457
00458 void H2_CollidRateRead( long int nColl )
00459 {
00460
00461
00462
00463
00464
00465
00466 const char* cdDATAFILE[N_X_COLLIDER] =
00467 {
00468 "H2_coll_H_07.dat",
00469 "H2_coll_He_LeBourlot.dat",
00470 "H2_coll_H2ortho.dat",
00471 "H2_coll_H2para.dat",
00472 "H2_coll_Hp.dat"
00473 };
00474
00475 FILE *ioDATA;
00476 char chLine[FILENAME_PATH_LENGTH_2];
00477 const char* chFilename;
00478 long int i, n1;
00479 long int iVibHi , iVibLo , iRotHi , iRotLo;
00480 long int magic_expect = -1;
00481
00482 DEBUG_ENTRY( "H2_CollidRateRead()" );
00483
00484 if( nColl == 0 )
00485 {
00486
00487 if( h2.lgH2_H_coll_07 )
00488 {
00489 magic_expect = 71106;
00490 chFilename = "H2_coll_H_07.dat";
00491 }
00492 else
00493 {
00494
00495 magic_expect = 20429;
00496 chFilename = "H2_coll_H_99.dat";
00497 }
00498 }
00499 else if( nColl == 1 )
00500 {
00501
00502 if( mole.lgH2_He_ORNL )
00503 {
00504
00505
00506 long int magic_found,
00507 magic_expect = 70513;
00508
00509
00510
00511
00512
00513
00514 chFilename = "H2_coll_He_ORNL.dat";
00515 if( (magic_found = H2_He_coll_init( chFilename )) != magic_expect )
00516 {
00517 fprintf(ioQQQ,"The H2 - He collision data file H2_coll_He_ORNL.dat does not have the correct magic number.\n");
00518 fprintf(ioQQQ,"I found %li but expected %li\n", magic_found , magic_expect );
00519
00520 cdEXIT(EXIT_FAILURE);
00521 }
00522 return;
00523 }
00524 else
00525 {
00526 magic_expect = 20429;
00527
00528 chFilename = "H2_coll_He_LeBourlot.dat";
00529 }
00530 }
00531
00532 else if( nColl == 2 )
00533 {
00534
00535 if( mole.lgH2_ORH2_ORNL )
00536 {
00537 long int magic_found,
00538 magic_expect = 80227;
00539
00540
00541
00542
00543
00544 char chPath[FILENAME_PATH_LENGTH_2];
00545 strcpy( chPath, "h2" );
00546 strcat( chPath, input.chDelimiter );
00547 strcat( chPath, "H2_coll_H2ortho_ORNL.dat" );
00548 if( (magic_found = H2_ORH2_coll_init( chPath )) != magic_expect )
00549 {
00550 fprintf(ioQQQ,"The H2 - H2 collision data file H2_coll_H2ortho_ORNL.dat does not have the correct magic number.\n");
00551 fprintf(ioQQQ,"I found %li but expected %li\n", magic_found , magic_expect );
00552
00553 cdEXIT(EXIT_FAILURE);
00554 }
00555 return;
00556 }
00557 else
00558 {
00559 magic_expect = 20429;
00560
00561 chFilename = "H2_coll_H2ortho_LeBourlot.dat";
00562 }
00563 }
00564
00565
00566 else if( nColl == 3 )
00567 {
00568
00569 if( mole.lgH2_PAH2_ORNL )
00570 {
00571 long int magic_found,
00572 magic_expect = 80227;
00573
00574 char chPath[FILENAME_PATH_LENGTH_2];
00575 strcpy( chPath, "h2" );
00576 strcat( chPath, input.chDelimiter );
00577 strcat( chPath, "H2_coll_H2para_ORNL.dat" );
00578 if( (magic_found = H2_PAH2_coll_init( chPath )) != magic_expect )
00579 {
00580 fprintf(ioQQQ,"The H2 - H2 collision data file H2_coll_H2para_ORNL.dat does not have the correct magic number.\n");
00581 fprintf(ioQQQ,"I found %li but expected %li\n", magic_found , magic_expect );
00582
00583 cdEXIT(EXIT_FAILURE);
00584 }
00585 return;
00586 }
00587 else
00588 {
00589 magic_expect = 20429;
00590
00591 chFilename = "H2_coll_H2para_LeBourlot.dat";
00592 }
00593 }
00594 else
00595 {
00596
00597 magic_expect = 20429;
00598 chFilename = cdDATAFILE[nColl];
00599 }
00600
00601 char chPath[FILENAME_PATH_LENGTH_2];
00602 strcpy( chPath, "h2" );
00603 strcat( chPath, input.chDelimiter );
00604 strcat( chPath, chFilename );
00605 ioDATA = open_data( chPath, "r" );
00606
00607
00608 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00609 {
00610 fprintf( ioQQQ, " H2_CollidRateRead could not read first line of %s\n", chFilename );
00611 cdEXIT(EXIT_FAILURE);
00612 }
00613
00614
00615 n1 = atoi( chLine );
00616
00617
00618
00619 if( n1 != magic_expect )
00620 {
00621 fprintf( ioQQQ,
00622 " H2_CollidRateRead: the version of %s is not the current version.\n", chFilename );
00623 fprintf( ioQQQ,
00624 " I expected to find the number %li and got %li instead.\n" ,
00625 magic_expect , n1 );
00626 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00627 cdEXIT(EXIT_FAILURE);
00628 }
00629
00630 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00631 {
00632 if( chLine[0]=='#' )
00633 continue;
00634 double a[3];
00635 sscanf(chLine,"%li\t%li\t%li\t%li\t%le\t%le\t%le",
00636 &iVibHi ,&iRotHi , &iVibLo , &iRotLo , &a[0],&a[1],&a[2] );
00637
00638
00639 ASSERT( iRotHi >= 0 && iVibHi >= 0 && iRotLo >= 0 && iVibLo >=0 );
00640
00641
00642
00643
00644 if( iVibHi > VIB_COLLID ||
00645 iVibLo > VIB_COLLID ||
00646 iRotHi > h2.nRot_hi[0][iVibHi] ||
00647 iRotLo > h2.nRot_hi[0][iVibLo])
00648 {
00649 iVibHi = -1;
00650 continue;
00651 }
00652
00653
00654 if( !( (iVibHi == iVibLo) && (iRotHi == iRotLo )) )
00655 {
00656
00657 ASSERT( (energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo] ) > 0. );
00658 for( i=0; i<3; ++i )
00659 {
00660 CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][i] = (realnum)a[i];
00661 }
00662
00663
00664
00665
00666 }
00667 }
00668 fclose( ioDATA );
00669
00670 return;
00671 }
00672
00673
00674
00675
00676
00677
00678 long int H2_He_coll_init(const char FILE_NAME_IN[])
00679 {
00680
00681 int i, j;
00682 long int magic;
00683 double par[N_H2_HE_FIT_PAR];
00684 char line[INPUT_LINE_LENGTH];
00685 int h2_i, h2_f, he_i, he_f;
00686 char quality, space;
00687
00688
00689 double error;
00690
00691 FILE *ifp;
00692
00693 DEBUG_ENTRY( "H2_He_coll_init()" );
00694
00695
00696
00697
00698 H2_He_coll_fit_par = (realnum ***)MALLOC(sizeof(realnum **)*(unsigned)nLevels_per_elec[0] );
00699 lgDefn_H2He_coll = (bool**)MALLOC(sizeof(bool *)*(unsigned)nLevels_per_elec[0] );
00700 for( i=1; i<nLevels_per_elec[0]; ++i )
00701 {
00702 lgDefn_H2He_coll[i] = (bool*)MALLOC(sizeof(bool)*(unsigned)i );
00703 H2_He_coll_fit_par[i] = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)i );
00704 for( j=0; j<i; ++j)
00705 H2_He_coll_fit_par[i][j] = (realnum *)MALLOC(sizeof(realnum)*(unsigned)N_H2_HE_FIT_PAR );
00706 }
00707
00708
00709 for( i=1; i<nLevels_per_elec[0]; i++ )
00710 {
00711 for( j=0; j<i; j++ )
00712 {
00713 lgDefn_H2He_coll[i][j] = false;
00714 }
00715 }
00716
00717
00718 char chPath[FILENAME_PATH_LENGTH_2];
00719 strcpy( chPath, "h2" );
00720 strcat( chPath, input.chDelimiter );
00721 strcat( chPath, FILE_NAME_IN );
00722 ifp = open_data( chPath, "r" );
00723
00724
00725 if( read_whole_line(line, (int)sizeof(line), ifp)==NULL )
00726 {
00727 printf("DISASTER H2_He_coll_init can't read first line of %s\n", FILE_NAME_IN);
00728 cdEXIT(EXIT_FAILURE);
00729 }
00730 sscanf(line, "%li", &magic);
00731
00732
00733 while( read_whole_line(line, (int)sizeof(line), ifp) != NULL )
00734 {
00735
00736 if( line[0]!='#' )
00737 {
00738 sscanf(line, "%i%i%i%i%c%c%c%c%lf%lf%lf%lf%lf%lf%lf%lf%lf", &h2_i,
00739 &h2_f, &he_i, &he_f,&space, &space, &space, &quality,
00740 &error, &par[0], &par[1], &par[2], &par[3], &par[4],
00741 &par[5], &par[6], &par[7]);
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763 if( par[3] < 0. )
00764 {
00765
00766 par[3] = -par[3];
00767 }
00768 if( par[6] < 0. )
00769 {
00770
00771 par[6] = -par[6];
00772 }
00773
00774
00775 --h2_f;
00776 --h2_i;
00777
00778 lgDefn_H2He_coll[h2_i][h2_f] = true;
00779 for( i=0; i<N_H2_HE_FIT_PAR; i++ )
00780
00781 H2_He_coll_fit_par[h2_i][h2_f][i] = (realnum)par[i];
00782 }
00783 }
00784 fclose(ifp);
00785 return magic;
00786 }
00787
00788
00789 long int H2_ORH2_coll_init(const char FILE_NAME_IN[] )
00790 {
00791
00792 int i, j;
00793 long int magic;
00794 double par[N_H2_ORH2_FIT_PAR];
00795 char line[INPUT_LINE_LENGTH];
00796 int h2_i, h2_f, h2ortho_i, h2ortho_f;
00797 char quality, space;
00798
00799
00800 double error;
00801
00802 FILE *ifp;
00803
00804 DEBUG_ENTRY( "H2_ORH2_coll_init()" );
00805
00806
00807
00808
00809 H2_ORH2_coll_fit_par = (realnum ***)MALLOC(sizeof(realnum **)*(unsigned)nLevels_per_elec[0] );
00810 lgDefn_H2ORH2_coll = (bool**)MALLOC(sizeof(bool *)*(unsigned)nLevels_per_elec[0] );
00811 for( i=1; i<nLevels_per_elec[0]; ++i )
00812 {
00813 lgDefn_H2ORH2_coll[i] = (bool*)MALLOC(sizeof(bool)*(unsigned)i );
00814 H2_ORH2_coll_fit_par[i] = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)i );
00815 for( j=0; j<i; ++j)
00816 H2_ORH2_coll_fit_par[i][j] = (realnum *)MALLOC(sizeof(realnum)*(unsigned)N_H2_ORH2_FIT_PAR );
00817 }
00818
00819
00820 for( i=1; i<nLevels_per_elec[0]; i++ )
00821 {
00822 for( j=0; j<i; j++ )
00823 {
00824 lgDefn_H2ORH2_coll[i][j] = false;
00825 }
00826 }
00827
00828
00829 ifp = open_data( FILE_NAME_IN, "r" );
00830
00831
00832 if( read_whole_line(line, (int)sizeof(line), ifp)==NULL )
00833 {
00834 printf("DISASTER H2_ORH2_coll_init can't read first line of %s\n", FILE_NAME_IN);
00835 cdEXIT(EXIT_FAILURE);
00836 }
00837 sscanf(line, "%li", &magic);
00838
00839
00840 while( read_whole_line(line, (int)sizeof(line), ifp) != NULL )
00841 {
00842
00843 if( line[0]!='#' )
00844 {
00845 sscanf(line, "%i%i%i%i%c%c%c%c%lf%lf%lf%lf%lf%lf%lf%lf%lf", &h2_i,
00846 &h2_f, &h2ortho_i, &h2ortho_f,&space, &space, &space, &quality,
00847 &error, &par[0], &par[1], &par[2], &par[3], &par[4],
00848 &par[5], &par[6], &par[7]);
00849
00850
00851 if( par[3] < 0. )
00852 {
00853
00854 par[3] = -par[3];
00855 }
00856 if( par[6] < 0. )
00857 {
00858
00859 par[6] = -par[6];
00860 }
00861
00862
00863 --h2_f;
00864 --h2_i;
00865
00866 lgDefn_H2ORH2_coll[h2_i][h2_f] = true;
00867 for( i=0; i<N_H2_ORH2_FIT_PAR; i++ )
00868
00869 H2_ORH2_coll_fit_par[h2_i][h2_f][i] = (realnum)par[i];
00870 }
00871 }
00872 fclose(ifp);
00873 return magic;
00874 }
00875
00876
00877 long int H2_PAH2_coll_init(const char FILE_NAME_IN[] )
00878 {
00879
00880 int i, j;
00881 long int magic;
00882 double par[N_H2_PAH2_FIT_PAR];
00883 char line[INPUT_LINE_LENGTH];
00884 int h2_i, h2_f, h2para_i, h2para_f;
00885 char quality, space;
00886
00887
00888 double error;
00889
00890 FILE *ifp;
00891
00892 DEBUG_ENTRY( "H2_PAH2_coll_init()" );
00893
00894
00895
00896
00897 H2_PAH2_coll_fit_par = (realnum ***)MALLOC(sizeof(realnum **)*(unsigned)nLevels_per_elec[0] );
00898 lgDefn_H2PAH2_coll = (bool**)MALLOC(sizeof(bool *)*(unsigned)nLevels_per_elec[0] );
00899 for( i=1; i<nLevels_per_elec[0]; ++i )
00900 {
00901 lgDefn_H2PAH2_coll[i] = (bool*)MALLOC(sizeof(bool)*(unsigned)i );
00902 H2_PAH2_coll_fit_par[i] = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)i );
00903 for( j=0; j<i; ++j)
00904 H2_PAH2_coll_fit_par[i][j] = (realnum *)MALLOC(sizeof(realnum)*(unsigned)N_H2_PAH2_FIT_PAR );
00905 }
00906
00907
00908 for( i=1; i<nLevels_per_elec[0]; i++ )
00909 {
00910 for( j=0; j<i; j++ )
00911 {
00912 lgDefn_H2PAH2_coll[i][j] = false;
00913 }
00914 }
00915
00916
00917 ifp = open_data( FILE_NAME_IN, "r" );
00918
00919 if( read_whole_line(line, (int)sizeof(line), ifp)==NULL )
00920 {
00921 printf("DISASTER H2_PAH2_coll_init can't read first line of %s\n", FILE_NAME_IN);
00922 cdEXIT(EXIT_FAILURE);
00923 }
00924 sscanf(line, "%li", &magic);
00925
00926
00927 while( read_whole_line(line, (int)sizeof(line), ifp) != NULL )
00928 {
00929
00930 if( line[0]!='#' )
00931 {
00932 sscanf(line, "%i%i%i%i%c%c%c%c%lf%lf%lf%lf%lf%lf%lf%lf%lf", &h2_i,
00933 &h2_f, &h2para_i, &h2para_f,&space, &space, &space, &quality,
00934 &error, &par[0], &par[1], &par[2], &par[3], &par[4],
00935 &par[5], &par[6], &par[7]);
00936
00937 if( par[3] < 0. )
00938 {
00939
00940 par[3] = -par[3];
00941 }
00942 if( par[6] < 0. )
00943 {
00944
00945 par[6] = -par[6];
00946 }
00947
00948
00949 --h2_f;
00950 --h2_i;
00951
00952 lgDefn_H2PAH2_coll[h2_i][h2_f] = true;
00953 for( i=0; i<N_H2_PAH2_FIT_PAR; i++ )
00954
00955 H2_PAH2_coll_fit_par[h2_i][h2_f][i] = (realnum)par[i];
00956 }
00957 }
00958 fclose(ifp);
00959 return magic;
00960 }
00961
00962
00963
00964
00965
00966 double H2_He_coll(int init, int final, double temp)
00967 {
00968
00969 double k, t3Plus1, b2, c2, f2, t3;
00970
00971 DEBUG_ENTRY( "H2_He_coll()" );
00972
00973
00974 if( temp<2 || temp > 1e8 )
00975 {
00976 k = -1;
00977 }
00978 else if( init <= final )
00979 {
00980 k = -1;
00981 }
00982
00983 else if( init < 0 || init >302 || final < 0 || final > 302 )
00984 {
00985 k = -1;
00986 }
00987
00988 else if( !lgDefn_H2He_coll[init][final] )
00989 {
00990 k = -1;
00991 }
00992
00993 else if( lgDefn_H2He_coll[init][final] )
00994 {
00995
00996
00997
00998
00999 double sum1 , sum2;
01000
01001
01002
01003
01005 temp = MIN2( 1e4 , temp );
01006
01007 t3 = temp/1e3;
01008
01009 t3Plus1 = t3+1;
01010 b2 = H2_He_coll_fit_par[init][final][1]/(H2_He_coll_fit_par[init][final][3]*t3+1);
01011 c2 = H2_He_coll_fit_par[init][final][2]/(t3Plus1*t3Plus1);
01012 {
01013 enum {DEBUG_LOC=false};
01014 if( DEBUG_LOC )
01015 {
01016 fprintf(ioQQQ,"bug H2 He coll\t%i %i %.3e %.3e %.3e \n",
01017 init,final,
01018 H2_He_coll_fit_par[init][final][5],
01019 H2_He_coll_fit_par[init][final][6]*t3+1,
01020 H2_He_coll_fit_par[init][final][7]
01021
01022 );
01023 }
01024 }
01025
01026 sum1 = H2_He_coll_fit_par[init][final][7] * log10( H2_He_coll_fit_par[init][final][6]*t3+1. );
01027
01028 if( fabs(sum1)< 38. )
01029 {
01030
01031
01032
01033
01034 f2 = H2_He_coll_fit_par[init][final][5]/pow(10. , sum1 );
01035 }
01036 else
01037 f2= 0.;
01038 {
01039 enum {DEBUG_LOC=false };
01040 if( DEBUG_LOC )
01041 {
01042 fprintf(ioQQQ,"bug H2 He coll\t%i %i %.3e %.3e %.3e %.3e %.3e sum %.3e %.3e \n",
01043 init,final,
01044 H2_He_coll_fit_par[init][final][0],
01045 b2,
01046 c2,
01047 H2_He_coll_fit_par[init][final][4],
01048 f2 ,
01049 H2_He_coll_fit_par[init][final][0]+b2+c2,
01050 H2_He_coll_fit_par[init][final][4]+f2);
01051 }
01052 }
01053
01054 sum1 = H2_He_coll_fit_par[init][final][0]+b2+c2;
01055 sum2 = H2_He_coll_fit_par[init][final][4]+f2;
01056 k = 0.;
01057 if( fabs(sum1) < 38. )
01058 k += pow(10., sum1 );
01059 if( fabs(sum2) < 38. )
01060 k += pow(10., sum2 );
01061
01062 if( k>1e-6 )
01063 {
01064
01065
01066 fprintf(ioQQQ,"PROBLEM H2-He rate coefficient hit cap, upper index of %i"
01067 " lower index of %i rate was %.2e resetting to 1e-6, Te=%e\n",
01068 init , final , k , phycon.te );
01069 k = 1e-6;
01070 }
01071 {
01072 enum {DEBUG_LOC=false};
01073 if( DEBUG_LOC )
01074 {
01075 fprintf(ioQQQ,"DEBUG H2 He rate %.3e \n",
01076 k);
01077 }
01078 }
01079 }
01080
01081 else
01082 {
01083 k = 99;
01084 }
01085 return k;
01086 }
01087
01088 double H2_ORH2_coll(int init, int final, double temp)
01089 {
01090
01091 double k, t, b2, c2, f2, t2;
01092
01093 DEBUG_ENTRY( "H2_ORH2_coll()" );
01094
01095
01096 if( temp<1 || temp > 1e4 )
01097 {
01098 k = -1;
01099 }
01100 else if( init <= final )
01101 {
01102 k = -1;
01103 }
01104
01105 else if( init < 0 || init >302 || final < 0 || final > 302 )
01106 {
01107 k = -1;
01108 }
01109
01110 else if( !lgDefn_H2ORH2_coll[init][final] )
01111 {
01112 k = -1;
01113 }
01114
01115 else if( lgDefn_H2ORH2_coll[init][final] )
01116 {
01117
01118
01119
01120
01121 double sum1 , sum2;
01122
01123
01124 t2 = temp/1e3;
01125
01126 t = t2+1;
01127 b2 = H2_ORH2_coll_fit_par[init][final][1]/(H2_ORH2_coll_fit_par[init][final][3]*t2+1);
01128 c2 = H2_ORH2_coll_fit_par[init][final][2]/(t*t);
01129 {
01130 enum {DEBUG_LOC=false};
01131 if( DEBUG_LOC )
01132 {
01133 fprintf(ioQQQ,"bug H2 H2ortho coll\t%i %i %.3e %.3e %.3e \n",
01134 init,final,
01135 H2_ORH2_coll_fit_par[init][final][5],
01136 H2_ORH2_coll_fit_par[init][final][6]*t2+1,
01137 H2_ORH2_coll_fit_par[init][final][7]
01138
01139 );
01140 }
01141 }
01142
01143 sum1 = H2_ORH2_coll_fit_par[init][final][7] * log10( H2_ORH2_coll_fit_par[init][final][6]*t2+1. );
01144
01145 if( fabs(sum1)< 38. )
01146 {
01147
01148 f2 = H2_ORH2_coll_fit_par[init][final][5]/pow(10. , sum1 );
01149 }
01150 else
01151 f2= 0.;
01152 {
01153 enum {DEBUG_LOC=false };
01154 if( DEBUG_LOC )
01155 {
01156 fprintf(ioQQQ,"bug H2 H2ortho coll\t%i %i %.3e %.3e %.3e %.3e %.3e sum %.3e %.3e \n",
01157 init,final,
01158 H2_ORH2_coll_fit_par[init][final][0],
01159 b2,
01160 c2,
01161 H2_ORH2_coll_fit_par[init][final][4],
01162 f2 ,
01163 H2_ORH2_coll_fit_par[init][final][0]+b2+c2,
01164 H2_ORH2_coll_fit_par[init][final][4]+f2);
01165 }
01166 }
01167
01168 sum1 = H2_ORH2_coll_fit_par[init][final][0]+b2+c2;
01169 sum2 = H2_ORH2_coll_fit_par[init][final][4]+f2;
01170 k = 0.;
01171 if( fabs(sum1) < 38. )
01172 k += pow(10., sum1 );
01173 if( fabs(sum2) < 38. )
01174 k += pow(10., sum2 );
01175
01176
01177 {
01178 enum {DEBUG_LOC=false};
01179 if( DEBUG_LOC )
01180 {
01181 fprintf(ioQQQ,"DEBUG H2 H2ortho rate %.3e \n",
01182 k);
01183 }
01184 }
01185 }
01186
01187 else
01188 {
01189 k = 99;
01190 }
01191 return k;
01192 }
01193
01194
01195 double H2_PAH2_coll(int init, int final, double temp)
01196 {
01197
01198 double k, t, b2, c2, f2, t2;
01199
01200 DEBUG_ENTRY( "H2_PAH2_coll()" );
01201
01202
01203 if( temp<1 || temp > 1e4 )
01204 {
01205 k = -1;
01206 }
01207 else if( init <= final )
01208 {
01209 k = -1;
01210 }
01211
01212 else if( init < 0 || init >302 || final < 0 || final > 302 )
01213 {
01214 k = -1;
01215 }
01216
01217 else if( !lgDefn_H2PAH2_coll[init][final] )
01218 {
01219 k = -1;
01220 }
01221
01222 else if( lgDefn_H2PAH2_coll[init][final] )
01223 {
01224
01225
01226
01227
01228 double sum1 , sum2;
01229
01230
01231 t2 = temp/1e3;
01232
01233 t = t2+1;
01234 b2 = H2_PAH2_coll_fit_par[init][final][1]/(H2_PAH2_coll_fit_par[init][final][3]*t2+1);
01235 c2 = H2_PAH2_coll_fit_par[init][final][2]/(t*t);
01236 {
01237 enum {DEBUG_LOC=false};
01238 if( DEBUG_LOC )
01239 {
01240 fprintf(ioQQQ,"bug H2 H2para coll\t%i %i %.3e %.3e %.3e \n",
01241 init,final,
01242 H2_PAH2_coll_fit_par[init][final][5],
01243 H2_PAH2_coll_fit_par[init][final][6]*t2+1,
01244 H2_PAH2_coll_fit_par[init][final][7]
01245
01246 );
01247 }
01248 }
01249
01250 sum1 = H2_PAH2_coll_fit_par[init][final][7] * log10( H2_PAH2_coll_fit_par[init][final][6]*t2+1. );
01251
01252 if( fabs(sum1)< 38. )
01253 {
01254
01255 f2 = H2_PAH2_coll_fit_par[init][final][5]/pow(10. , sum1 );
01256 }
01257 else
01258 f2= 0.;
01259 {
01260 enum {DEBUG_LOC=false };
01261 if( DEBUG_LOC )
01262 {
01263 fprintf(ioQQQ,"bug H2 H2para coll\t%i %i %.3e %.3e %.3e %.3e %.3e sum %.3e %.3e \n",
01264 init,final,
01265 H2_PAH2_coll_fit_par[init][final][0],
01266 b2,
01267 c2,
01268 H2_PAH2_coll_fit_par[init][final][4],
01269 f2 ,
01270 H2_PAH2_coll_fit_par[init][final][0]+b2+c2,
01271 H2_PAH2_coll_fit_par[init][final][4]+f2);
01272 }
01273 }
01274
01275 sum1 = H2_PAH2_coll_fit_par[init][final][0]+b2+c2;
01276 sum2 = H2_PAH2_coll_fit_par[init][final][4]+f2;
01277 k = 0.;
01278 if( fabs(sum1) < 38. )
01279 k += pow(10., sum1 );
01280 if( fabs(sum2) < 38. )
01281 k += pow(10., sum2 );
01282
01283
01284 {
01285 enum {DEBUG_LOC=false};
01286 if( DEBUG_LOC )
01287 {
01288 fprintf(ioQQQ,"DEBUG H2 H2para rate %.3e \n",
01289 k);
01290 }
01291 }
01292 }
01293
01294 else
01295 {
01296 k = 99;
01297 }
01298 return k;
01299 }
01300