00001
00002
00003
00004 #include "cddefines.h"
00005 #include "ionbal.h"
00006 #include "dense.h"
00007 #include "phycon.h"
00008 #include "save.h"
00009 #include "atmdat.h"
00010
00011
00012
00013
00014
00015
00016
00017
00018 #ifdef USENEW
00019 STATIC double dr_suppress(
00020
00021
00022 long int atomic_number,
00023
00024 long int ionic_charge,
00025
00026 double eden,
00027
00028 double T )
00029 {
00030
00031
00032
00033 const double A = 12.4479;
00034 const double mu = 0.46665;
00035 const double w = 4.96916;
00036 const double y_c0 = 0.5498;
00037
00038
00039
00040
00041 const double T_0 = 1.e5;
00042 const double q_0 = 3.;
00043
00044
00045
00046
00047
00048 const double c = 3.0;
00049
00050 double s, snew, y_c, E_c, E_c0, x, g_x;
00051 long int iso_sequence;
00052
00053 eden = log10(eden);
00054 iso_sequence = atomic_number - ionic_charge;
00055
00056 y_c = y_c0 + log10( pow( (ionic_charge/q_0), 7. ) * sqrt( T/T_0 ) );
00057
00058
00059 if( eden >= y_c )
00060 {
00061 s = (A/(PI*w)) * ( mu/( 1. + pow((eden-y_c)/w, 2.) ) +
00062 (1. - mu) * sqrt(PI*LN_TWO) * exp( -LN_TWO *
00063 pow((eden-y_c)/w, 2.) ) );
00064 }
00065 else
00066 {
00067 s = (A/(PI*w)) * ( mu + (1.- mu) * sqrt(PI*LN_TWO) );
00068 }
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079 if( iso_sequence == 3 )
00080 {
00081 E_c = 2.08338 + 19.1356*(ionic_charge/10.) + 0.974 *
00082 pow( ionic_charge/10., 2. ) - 0.309032*pow( ionic_charge/10., 3. ) +
00083 0.419951*pow( ionic_charge/10., 4. );
00084 }
00085 else if( iso_sequence == 4 )
00086 {
00087 E_c = 5.56828 + 34.6774*(ionic_charge/10.) + 1.005 *
00088 pow( ionic_charge/10., 2. ) - 0.994177*pow( ionic_charge/10., 3. ) +
00089 0.726053*pow( ionic_charge/10., 4. );
00090 }
00091 else if( iso_sequence == 7 )
00092 {
00093 E_c = 10.88361 + 39.7851*(ionic_charge/10.) + 0.423 *
00094 pow( ionic_charge/10., 2. ) - 0.310368*pow( ionic_charge/10., 3. ) +
00095 0.937186*pow( ionic_charge/10., 4. );
00096 }
00097 else if( iso_sequence == 11 )
00098 {
00099 E_c = 2.17262 + 22.5038*(ionic_charge/10.) - 1.227*pow( ionic_charge/10., 2. ) +
00100 0.801291*pow( ionic_charge/10., 3. ) +
00101 0.0434168*pow( ionic_charge/10., 4. );
00102 }
00103 else if( iso_sequence == 1 || iso_sequence == 2 || iso_sequence == 10 )
00104 {
00105
00106
00107 E_c = 1.e10;
00108 }
00109 else
00110 {
00111
00112 E_c = 0.0;
00113
00114
00115 }
00116
00117
00118 E_c0 = 2.08338 + 19.1356*(q_0/10.) + 0.974 *
00119 pow( (q_0/10.), 2. ) - 0.309032 *
00120 pow( (q_0/10.), 3. ) + 0.419951 *
00121 pow( (q_0/10.), 4. );
00122
00123
00124 x = ( (E_c*EVDEGK)/(c*T) - (E_c0*EVDEGK)/(c*T_0) );
00125
00126 if( x > 1 )
00127 {
00128 g_x = x;
00129 }
00130 else if( x >= 0 && x <= 1 )
00131 {
00132 g_x = x*x;
00133 }
00134 else
00135 {
00136 g_x = 0.0;
00137 }
00138
00139
00140
00141
00142 snew = 1. + (s-1.)*exp(-g_x);
00143
00144 return snew;
00145 ASSERT( snew >=0. && snew <= 1. );
00146 }
00147 #endif
00148
00149 void atmdat_DielSupres(void)
00150 {
00151 long int i;
00152
00153 DEBUG_ENTRY( "atmdat_DielSupres()" );
00154
00155
00156 if( ionbal.lgSupDie[0] )
00157 {
00158 for( i=0; i < LIMELM; i++ )
00159 {
00160 # ifdef USENEW
00161 ionbal.DielSupprs[0][i] = (realnum)dr_suppress( i+1, 3, dense.eden, phycon.te );
00162 # else
00163
00164
00165 double effden = dense.eden/(phycon.sqrte/122.47);
00166
00167
00168 effden /= powi((realnum)(i+1)/3.,7);
00169
00170 ionbal.DielSupprs[0][i] = (realnum)(1.-0.092*log10(effden));
00171 ionbal.DielSupprs[0][i] = (realnum)MIN2(1.,ionbal.DielSupprs[0][i]);
00172 ionbal.DielSupprs[0][i] = (realnum)MAX2(0.08,ionbal.DielSupprs[0][i]);
00173 # endif
00174 }
00175 }
00176
00177 else
00178 {
00179 for( i=0; i < LIMELM; i++ )
00180 {
00181 ionbal.DielSupprs[0][i] = 1.;
00182 }
00183 }
00184
00185
00186
00187 if( ionbal.lgSupDie[1] )
00188 {
00189 for( i=0; i < LIMELM; i++ )
00190 {
00191
00192 ionbal.DielSupprs[1][i] = ionbal.DielSupprs[0][i];
00193 }
00194 }
00195 else
00196 {
00197 for( i=0; i < LIMELM; i++ )
00198 {
00199 ionbal.DielSupprs[1][i] = 1.;
00200 }
00201 }
00202
00203
00204
00205 if( save.lgioRecom )
00206 {
00207 fprintf( save.ioRecom, " atmdat_DielSupres finds following dielectronic"
00208 " recom suppression factors.\n" );
00209 fprintf( save.ioRecom, " Z fac \n" );
00210 for( i=0; i < LIMELM; i++ )
00211 {
00212 fprintf( save.ioRecom, "%3ld %10.3e %10.3e\n", i+1,
00213 ionbal.DielSupprs[0][i], ionbal.DielSupprs[1][i] );
00214 }
00215 fprintf( save.ioRecom, "\n");
00216 }
00217 return;
00218 }