00001
00002
00003
00004
00005 #include "cddefines.h"
00006 #include "atmdat.h"
00007 #include "phycon.h"
00008 #include "dense.h"
00009 #include "taulines.h"
00010 #include "input.h"
00011 #include "h2.h"
00012 #include "h2_priv.h"
00013
00014 STATIC realnum GbarRateCoeff( long nColl, double ediff );
00015
00016
00017 void diatomics::H2_CollidRateEvalAll( void )
00018 {
00019 DEBUG_ENTRY( "H2_CollidRateEvalAll()" );
00020
00021 long int numb_coll_trans = 0;
00022
00023 if( nTRACE >= n_trace_full )
00024 fprintf(ioQQQ,"%s set collision rates\n", label.c_str());
00025
00026
00027
00028
00029
00030 H2_coll_dissoc_rate_coef[0][0] = 0.;
00031 H2_coll_dissoc_rate_coef_H2[0][0] = 0.;
00032 for( long ipHi=0; ipHi<nLevels_per_elec[0]; ++ipHi )
00033 {
00034
00035 long iVibHi = ipVib_H2_energy_sort[ipHi];
00036 long iRotHi = ipRot_H2_energy_sort[ipHi];
00037
00038
00039
00040
00041 double energy = H2_DissocEnergies[0] - states[ipHi].energy().WN();
00042 ASSERT( energy > 0. );
00043
00044 H2_coll_dissoc_rate_coef[iVibHi][iRotHi] =
00045 1e-14f * (realnum)sexp(energy/phycon.te_wn) * lgColl_dissoc_coll;
00046
00047
00048
00049 H2_coll_dissoc_rate_coef_H2[iVibHi][iRotHi] =
00050 1e-11f * (realnum)sexp(energy/phycon.te_wn) * lgColl_dissoc_coll;
00051
00052
00053
00054
00055
00056
00057 for( long ipLo=0; ipLo<ipHi; ++ipLo )
00058 {
00059 long iVibLo = ipVib_H2_energy_sort[ipLo];
00060 long iRotLo = ipRot_H2_energy_sort[ipLo];
00061
00062 ASSERT( states[ipHi].energy().WN() - states[ipLo].energy().WN() > 0. );
00063
00064
00065 ++numb_coll_trans;
00066 for( long nColl=0; nColl<N_X_COLLIDER; ++nColl )
00067 {
00068 realnum newrate = H2_CollidRateEvalOne( iVibHi,iRotHi,iVibLo,iRotLo,
00069 ipHi , ipLo , nColl, phycon.te );
00070 CollRateCoeff[ipHi][ipLo][nColl] = newrate;
00071 }
00072 }
00073 }
00074
00075 fixit();
00076 #if 0
00077
00078
00079
00080
00081
00082
00083 if( phycon.te < 100. )
00084 {
00085
00086
00087 H2_CollRate[0][1][0][0][0] *= (realnum)(exp(-(3900-170.5)*(1./phycon.te - 1./100.)));
00088 H2_CollRate[0][3][0][0][0] *= (realnum)(exp(-(3900-1015.1)*(1./phycon.te - 1./100.)));
00089 H2_CollRate[0][2][0][1][0] *= (realnum)(exp(-(3900-339.3)*(1./phycon.te - 1./100.)));
00090 }
00091 #endif
00092
00093 if( nTRACE >= n_trace_full )
00094 fprintf(ioQQQ,
00095 " collision rates updated for new temp, number of trans is %li\n",
00096 numb_coll_trans);
00097
00098 return;
00099 }
00100
00101
00102 realnum diatomics::H2_CollidRateEvalOne(
00103
00104
00105 long iVibHi, long iRotHi,long iVibLo,
00106
00107 long iRotLo, long ipHi , long ipLo ,
00108
00109 long nColl,
00110 double te_K )
00111 {
00112 DEBUG_ENTRY( "H2_CollidRateEvalOne()" );
00113
00114
00115
00116
00117 realnum rate = InterpCollRate( RateCoefTable[nColl], ipHi, ipLo, te_K );
00118
00119
00120
00121
00122 if( (rate==0.f) && lgColl_gbar &&
00123 ( nColl==4 ||
00124 ( ( nColl==0 || nColl==1 ) &&
00125 H2_lgOrtho[0][iVibHi][iRotHi]-H2_lgOrtho[0][iVibLo][iRotLo]==0 ) ) )
00126 {
00127 double ediff = states[ipHi].energy().WN() - states[ipLo].energy().WN();
00128 rate = GbarRateCoeff( nColl, ediff );
00129 }
00130
00131 rate *= lgColl_deexec_Calc;
00132
00133
00134 if( !lgH2_ortho_para_coll_on &&
00135 (H2_lgOrtho[0][iVibHi][iRotHi]-H2_lgOrtho[0][iVibLo][iRotLo]) )
00136 rate = 0.;
00137
00138 if( lgH2_NOISE )
00139 rate *= CollRateErrFac[ipHi][ipLo][nColl];
00140
00141 return rate;
00142 }
00143
00144 STATIC realnum GbarRateCoeff( long nColl, double ediff )
00145 {
00146
00147
00148 double gbarcoll[N_X_COLLIDER][3] =
00149 {
00150 {-9.9265 , -0.1048 , 0.456 },
00151 {-8.281 , -0.1303 , 0.4931 },
00152 {-10.0357, -0.0243 , 0.67 },
00153 {-8.6213 , -0.1004 , 0.5291 },
00154 {-9.2719 , -0.0001 , 1.0391 }
00155 };
00156
00157
00158
00159 ediff = MAX2(100., ediff );
00160
00161
00162
00163 realnum rate = (realnum)pow(10. ,
00164 gbarcoll[nColl][0] + gbarcoll[nColl][1] *
00165 pow(ediff,gbarcoll[nColl][2]) );
00166
00167 return rate;
00168 }
00169
00170 void diatomics::H2_CollidRateRead( long int nColl )
00171 {
00172 DEBUG_ENTRY( "H2_CollidRateRead()" );
00173
00174 if( coll_source[nColl].filename.size() == 0 && coll_source[nColl].magic == 0 )
00175 return;
00176
00177 const char* chFilename = coll_source[nColl].filename.c_str();
00178 long magic_expect = coll_source[nColl].magic;
00179 char chLine[INPUT_LINE_LENGTH];
00180 char chPath[FILENAME_PATH_LENGTH_2];
00181 strcpy( chPath, path.c_str() );
00182 strcat( chPath, input.chDelimiter );
00183 strcat( chPath, chFilename );
00184 FILE *ioDATA = open_data( chPath, "r" );
00185
00186
00187 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00188 {
00189 fprintf( ioQQQ, " H2_CollidRateRead could not read first line of %s\n", chFilename );
00190 cdEXIT(EXIT_FAILURE);
00191 }
00192
00193
00194 long n1 = atoi( chLine );
00195 if( n1 != magic_expect )
00196 {
00197 fprintf( ioQQQ, " H2_CollidRateRead: the version of %s is not the current version.\n", chFilename );
00198 fprintf( ioQQQ, " I expected to find the number %li and got %li instead.\n", magic_expect, n1 );
00199 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00200 cdEXIT(EXIT_FAILURE);
00201 }
00202
00203 FunctPtr func = new FunctDiatoms( *this );
00204 ReadCollisionRateTable( RateCoefTable[nColl], ioDATA, func, nLevels_per_elec[0] );
00205 delete func;
00206
00207 fclose( ioDATA );
00208
00209 return;
00210 }
00211
00212 void diatomics::GetIndices( long& ipHi, long& ipLo, const char* chLine, long& i ) const
00213 {
00214 bool lgEOL;
00215 long iVibHi = (long)FFmtRead( chLine, &i, strlen(chLine), &lgEOL );
00216 long iRotHi = (long)FFmtRead( chLine, &i, strlen(chLine), &lgEOL );
00217 long iVibLo = (long)FFmtRead( chLine, &i, strlen(chLine), &lgEOL );
00218 long iRotLo = (long)FFmtRead( chLine, &i, strlen(chLine), &lgEOL );
00219 ASSERT( iRotHi >= 0 && iVibHi >= 0 && iRotLo >= 0 && iVibLo >=0 );
00220
00221
00222
00223
00224 if( ( iVibHi > nVib_hi[0] || iVibLo > nVib_hi[0] ) ||
00225 ( iRotHi < Jlowest[0] || iRotLo < Jlowest[0] ) ||
00226 ( iRotHi > nRot_hi[0][iVibHi] || iRotLo > nRot_hi[0][iVibLo] ) ||
00227
00228 ( iVibHi == iVibLo && iRotHi == iRotLo ) )
00229 {
00230 ipHi = -1;
00231 ipLo = -1;
00232 return;
00233 }
00234
00235 ipHi = ipEnergySort[0][iVibHi][iRotHi];
00236 ipLo = ipEnergySort[0][iVibLo][iRotLo];
00237 if( ipHi < ipLo )
00238 swap( ipHi, ipLo );
00239
00240 return;
00241 }
00242