00001 /* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 #include "cddefines.h" 00004 #include "deuterium.h" 00005 #include "dense.h" 00006 #include "mole.h" 00007 00008 t_deuterium deut; 00009 00010 void ScaleDensitiesDeuterium( const realnum &factor ) 00011 { 00012 deut.gas_phase *= factor; 00013 deut.xMolecules *= factor; 00014 deut.xIonDense[0] *= factor; 00015 deut.xIonDense[1] *= factor; 00016 return; 00017 } 00018 00019 void InitDeuteriumIonization() 00020 { 00021 deut.xIonDense[0] = deut.gas_phase; 00022 deut.xIonDense[1] = 0.; 00023 return; 00024 } 00025 00026 void SetDeuteriumIonization( const double &xNeutral, const double &xIonized ) 00027 { 00028 if( !deut.lgElmtOn ) 00029 return; 00030 00031 total_molecule_deut( deut.xMolecules ); 00032 00033 realnum total = deut.gas_phase - deut.xMolecules; 00034 00035 fixit(); // try to let D ionization be independent of H 00036 #if 1 00037 realnum neut = total * xNeutral/( xNeutral + xIonized ); 00038 realnum ionz = total * xIonized/( xNeutral + xIonized ); 00039 #else 00040 realnum src = mole.findrk("D+,e-=>D,PHOTON") * dense.eden; 00041 realnum snk = mole.findrk("D,PHOTON=>D+,e-"); 00042 realnum ion_ratio = xIonized/xNeutral; 00043 if( src > SMALLFLOAT ) 00044 ion_ratio = snk/src; 00045 realnum neut = total / ( ion_ratio + 1.f ); 00046 realnum ionz = total * ion_ratio / ( ion_ratio + 1.f ); 00047 #endif 00048 if( total > 1e-4 * deut.gas_phase ) 00049 { 00050 deut.xIonDense[0] = neut; 00051 deut.xIonDense[1] = ionz; 00052 } 00053 00054 //if( iteration >= 3 ) 00055 // fprintf( ioQQQ, "DEBUGGG SetDeuteriumIonization %li\t%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%.12e\t%.12e\n", 00056 // nzone, deut.gas_phase, deut.xMolecules, total, old0, deut.xIonDense[0], old1, deut.xIonDense[1] ); 00057 00058 return; 00059 } 00060 00061 void SetDeuteriumFractionation( const realnum &frac ) 00062 { 00063 if( !deut.lgElmtOn ) 00064 return; 00065 deut.fractionation = frac; 00066 return; 00067 } 00068 00069 void SetGasPhaseDeuterium( const realnum &Hdensity ) 00070 { 00071 if( !deut.lgElmtOn ) 00072 return; 00073 deut.gas_phase = Hdensity * deut.fractionation; 00074 } 00075