cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*eden_sum sum free electron density over all species, sets variable erredn.EdenTrue 00004 *called by ConvBase - ConvEdenIoniz actually updates the electron density 00005 * returns 0 if all is ok, 1 if need to abort calc */ 00006 #include "cddefines.h" 00007 #include "hmi.h" 00008 #include "trace.h" 00009 #include "grainvar.h" 00010 #include "rfield.h" 00011 #include "mole.h" 00012 #include "dense.h" 00013 #include "conv.h" 00014 00015 /*eden_sum sum free electron density over all species, sets variable erredn.EdenTrue 00016 *called by ConvBase - ConvEdenIoniz actually updates the electron density 00017 * returns 0 if all is ok, 1 if need to abort calc */ 00018 int eden_sum(void) 00019 { 00020 long int i, 00021 ion, 00022 nelem; 00023 00024 realnum sum_all_ions , 00025 sum_metals , 00026 hmole_eden, 00027 eden_ions[LIMELM]; 00028 00029 DEBUG_ENTRY( "eden_sum()" ); 00030 00031 /* EdenExtra is normally zero, set with EDEN command, to add extra e- */ 00032 dense.EdenTrue = dense.EdenExtra; 00033 00034 /* sum over all ions */ 00035 sum_all_ions = 0.; 00036 sum_metals = 0.; 00037 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ ) 00038 { 00039 if( nelem==ipLITHIUM ) 00040 sum_metals = 0.; 00041 eden_ions[nelem] = dense.xIonDense[nelem][1]; 00042 /* >>chng 02 jul 15, upper limit now includes H-like, so all species 00043 * in this one loop (except hydrogen and helium ) */ 00044 /* >>chng 02 apr 26, upper limit had been nelem, not nelem+1, so H-like 00045 * metals were missed - ok for H and He since not part of this loop before */ 00046 for( ion=2; ion <= nelem+1; ion++ ) 00047 { 00048 /* >>chng 96 oct 27, save all contributors to electron density */ 00049 eden_ions[nelem] += ion*dense.xIonDense[nelem][ion]; 00050 } 00051 00052 sum_all_ions += eden_ions[nelem]; 00053 sum_metals += eden_ions[nelem]; 00054 } 00055 dense.EdenTrue += sum_all_ions; 00056 ASSERT( dense.EdenTrue >= 0. ); 00057 00058 /* add on molecules */ 00059 co.comole_eden = 0.; 00060 for( i=0; i < mole.num_comole_calc; i++ ) 00061 { 00062 if(COmole[i]->n_nuclei != 1) 00063 { 00064 co.comole_eden += COmole[i]->hevmol*COmole[i]->nElec; 00065 } 00066 } 00067 00068 /* electrons contributed by molecules */ 00069 dense.EdenTrue += co.comole_eden; 00070 ASSERT( dense.EdenTrue >= 0. ); 00071 00072 /* >>chng 03 nov 28, loop over all H molecules, had just been H- */ 00073 hmole_eden = 0.; 00074 for(i=0;i<N_H_MOLEC;i++) 00075 { 00076 /* hmi.nElectron is zero for H+ since counted as an ion, -1 for H-, etc */ 00077 hmole_eden += hmi.Hmolec[i]*hmi.nElectron[i]; 00078 } 00079 /* >>chng 01 jan 08, following logic added to stop H- from forcing 00080 * negative electron density during very first search, when H- is badly off */ 00081 /* >>chng 03 nov 28, add test for lgSearch, had been absent */ 00082 /* >>chng 04 feb 20, recoil_ion.in shows this important, 00083 * update test so prevent neg eden */ 00084 if( (-hmole_eden) > dense.EdenTrue/4. && conv.lgSearch ) 00085 { 00086 /*dense.EdenTrue += MIN2(dense.EdenTrue*0.05,hmole_eden);*/ 00087 dense.EdenTrue *= 0.9; 00088 } 00089 else if( (-hmole_eden) > dense.EdenTrue ) 00090 { 00091 /* >>chng 04 mar 14, add this branch to prevent 00092 * negative electron density. occurred in pdr 00093 * with large h2, after hmole failure */ 00094 fprintf(ioQQQ," PROBLEM sum eden from hmole too neg, set to limt. EdenTrue:%.3e hmole_eden:%.3e \n", 00095 dense.EdenTrue, 00096 hmole_eden); 00097 dense.EdenTrue = dense.EdenTrue/2.; 00098 } 00099 else 00100 { 00101 /* account for electrons on all H-bearing molecules */ 00102 dense.EdenTrue += hmole_eden; 00103 } 00104 ASSERT(dense.EdenTrue >= 0. ); 00105 00106 /* this variable is set with the set eden command, 00107 * is supposed to override physical electron density */ 00108 if( dense.EdenSet > 0. ) 00109 { 00110 dense.EdenTrue = dense.EdenSet; 00111 dense.eden_from_metals = 1.; 00112 } 00113 else 00114 { 00115 /* fraction of electrons from si, s, mg, al */ 00116 /*dense.eden_from_metals = sum_all_ions / SDIV(dense.EdenTrue);*/ 00117 dense.eden_from_metals = sum_metals / SDIV(dense.EdenTrue); 00118 /*fprintf(ioQQQ," debuggg ionfrac %.2f %.3f\n",fnzone,dense.eden_from_metals );*/ 00119 } 00120 /* dense.eden itself is actually changed in ConvEdenIoniz */ 00121 00122 /* >>chng 00 dec 19, include electrons on grains in total sum */ 00123 /* >>chng 02 jul 15, only add grain electrons when we have stable soln - 00124 * everett.in had very negative grain elec total, so drove total eden negative, 00125 * but due to bad inital electron density - do not add grain electrons until we 00126 * are close to a solution */ 00127 /* >>chng 04 sep 26, allow nEdenTrue to become -ve */ 00128 if( dense.EdenTrue+gv.TotalEden*gv.lgGrainElectrons >= 0. ) 00129 { 00130 /* gv.lgGrainElectrons - should grain electron source/sink be included in overall electron sum? 00131 * default is true, set false with no grain electrons command */ 00132 dense.EdenTrue += gv.TotalEden*gv.lgGrainElectrons; 00133 } 00134 else 00135 { 00136 /* >>chng 05 jan 24, only print if not in search phase */ 00137 if( !conv.lgSearch ) 00138 fprintf(ioQQQ, 00139 " PROBLEM eden grain neg limt %.2f eden %.4e new eden bef grn %.4e" 00140 "grain eden %.4e set edentrue to %.4e Search?%c\n", 00141 fnzone, 00142 dense.eden, 00143 dense.EdenTrue, 00144 gv.TotalEden, 00145 dense.eden*0.9, 00146 TorF(conv.lgSearch)); 00147 00148 /*dense.EdenTrue = dense.eden*0.9;*/ 00149 dense.EdenTrue += gv.TotalEden*gv.lgGrainElectrons; 00150 } 00151 00152 if( trace.lgTrace || trace.lgESOURCE ) 00153 { 00154 fprintf( ioQQQ, 00155 " eden_sum zn:%.2f current:%.4e new true:%.4e ions:%.4e comole:%.4e hmole:%.4e grain:%.4e extra:%.4e sum:%.4e LaOTS:%.4e\n", 00156 fnzone , 00157 dense.eden , 00158 dense.EdenTrue , 00159 sum_all_ions , 00160 co.comole_eden , 00161 hmole_eden , 00162 gv.TotalEden*gv.lgGrainElectrons, 00163 dense.EdenExtra , 00164 sum_all_ions + co.comole_eden + hmole_eden + gv.TotalEden*gv.lgGrainElectrons + dense.EdenExtra, 00165 rfield.otslin[2182] ); 00166 00167 if( trace.lgNeBug ) 00168 { 00169 for(nelem=ipHYDROGEN; nelem < LIMELM; nelem++) 00170 { 00171 if( nelem==0 ) 00172 { 00173 fprintf( ioQQQ, " eden_sum H -Ne:" ); 00174 } 00175 else if( nelem==10 ) 00176 { 00177 fprintf( ioQQQ, "\n eden_sum Na-Ca:" ); 00178 } 00179 else if( nelem==20 ) 00180 { 00181 fprintf( ioQQQ, "\n eden_sum Sc-Zn:" ); 00182 } 00183 fprintf( ioQQQ, " %.4e", eden_ions[nelem] ); 00184 if( nelem==29 ) 00185 { 00186 fprintf( ioQQQ, "\n" ); 00187 } 00188 /*if( nelem==9 || nelem==19 || nelem==29 ) 00189 fprintf( ioQQQ, "\n " );*/ 00190 } 00191 } 00192 } 00193 00194 /* abort if negative electron density */ 00195 /*>>chng 04 sep 25, allow negative true elec density so solver can work */ 00196 if( 0 && dense.EdenTrue < 0. ) 00197 { 00198 fprintf( ioQQQ, "eden_sum finds non-positive electron density.\n" ); 00199 fprintf( ioQQQ, 00200 " eden_sum: EdenTrue to%12.4e fm%12.4e Ne(H):%10.2e Ne(He):%10.2e Ne(C):\n", 00201 dense.EdenTrue, 00202 dense.eden, 00203 dense.xIonDense[ipHYDROGEN][1], 00204 dense.xIonDense[ipHELIUM][1] + 2.*dense.xIonDense[ipHELIUM][2] ); 00205 ShowMe(); 00206 cdEXIT(EXIT_FAILURE); 00207 } 00208 else if( dense.EdenTrue == 0. ) 00209 { 00210 fprintf( ioQQQ, "\nDISASTER PROBLEM eden_sum finds an electron density of zero, this is unphysical.\n" ); 00211 fprintf( ioQQQ, "There appears to be no source of ionization.\n"); 00212 fprintf( ioQQQ, "Consider adding background cosmic rays with COSMIC RAY BACKGROUND and BACKGROUND commands.\n\n"); 00213 lgAbort = true; 00214 00215 return 1; 00216 } 00217 00218 { 00219 /*@-redef@*/ 00220 enum {DEBUG_LOC=false}; 00221 /*@+redef@*/ 00222 if( DEBUG_LOC ) 00223 { 00224 fprintf(ioQQQ,"esumdebugg\t%li\t%.2e\t%.2e\n", 00225 nzone, 00226 dense.eden, dense.EdenTrue); 00227 } 00228 } 00229 00230 /* >>chng 05 jan 05, don't let elec den be zero - logs are taken */ 00231 dense.eden = MAX2( SMALLFLOAT , dense.eden ); 00232 00233 return 0; 00234 }