Changeset 2086 for branches/newmole
- Timestamp:
- 05/19/08 12:39:03 (8 months ago)
- Location:
- branches/newmole/source
- Files:
-
- 5 modified
-
iso_create.cpp (modified) (5 diffs)
-
iso_level.cpp (modified) (2 diffs)
-
newton_step.cpp (modified) (2 diffs)
-
parse_caseb.cpp (modified) (2 diffs)
-
radius_increment.cpp (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/newmole/source/iso_create.cpp
r2023 r2086 331 331 StatesElem[ipISO][nelem][ipHi].lifetime = SMALLFLOAT; 332 332 333 long ipLoStart = 0; 334 if( opac.lgCaseB && L_(ipHi)==1 && (ipISO==ipH_LIKE || S_(ipHi)==1) ) 335 ipLoStart = 1; 336 337 for( ipLo=ipLoStart; ipLo < ipHi; ipLo++ ) 333 for( ipLo=0; ipLo < ipHi; ipLo++ ) 338 334 { 339 335 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA ) … … 1272 1268 for( ipHi=1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ ) 1273 1269 { 1270 double SumAs = 0.; 1271 1274 1272 /** Cascade probabilities are as defined in Robbins 68, 1275 1273 * generalized here for cascade probability for any iso sequence. … … 1294 1292 for( ipLo=ipLoStart; ipLo<ipHi; ipLo++ ) 1295 1293 { 1294 SumAs += Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul; 1295 } 1296 1297 for( ipLo=ipLoStart; ipLo<ipHi; ipLo++ ) 1298 { 1296 1299 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA ) 1297 1300 { … … 1303 1306 iso.CascadeProb[ipISO][nelem][ipHi][ipLo] = 0.; 1304 1307 iso.BranchRatio[ipISO][nelem][ipHi][ipLo] = 1305 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul * 1306 StatesElem[ipISO][nelem][ipHi].lifetime; 1308 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul / SumAs; 1307 1309 1308 1310 ASSERT( iso.BranchRatio[ipISO][nelem][ipHi][ipLo] <= 1.0000001 ); … … 1965 1967 StatesElem[ipISO][nelem][ipHi].lifetime = SMALLFLOAT; 1966 1968 1967 long ipLoStart=0; 1968 if( opac.lgCaseB && L_(ipHi)==1 && (ipISO==ipH_LIKE || S_(ipHi)==1) ) 1969 ipLoStart=1; 1970 1971 for( long ipLo=ipLoStart; ipLo < ipHi; ipLo++ ) 1969 for( long ipLo=0; ipLo < ipHi; ipLo++ ) 1972 1970 { 1973 1971 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA ) -
branches/newmole/source/iso_level.cpp
r2080 r2086 334 334 } 335 335 336 for( long ion=0; ion<=nelem+1; ++ion ) 337 { 338 if( ion!=nelem-ipISO ) 339 { 340 /* recombination must be multiplied by a ratio of densities to get proper rate. */ 341 source += mole.xMoleChTrRate[nelem][ion][nelem-ipISO] * 342 dense.xIonDense[nelem][ion] ; 343 /* take ionization out of every level. */ 344 sink += mole.xMoleChTrRate[nelem][nelem-ipISO][ion]; 345 } 346 } 347 336 348 /* add in source and sink terms from molecular network. */ 337 if( nelem-ipISO < N_MOLE_ION)349 if( conv.nTotalIoniz > 1 || iteration > 1) 338 350 { 339 351 source += mole.source[nelem][nelem-ipISO]; … … 350 362 } 351 363 } 364 } 365 366 if( nelem-ipISO >= 1 && ionbal.RateIonizTot[nelem][nelem-ipISO-1] > 0.) 367 { 368 /* now add ionization from next lower stage */ 369 source += ionbal.RateIonizTot[nelem][nelem-ipISO-1]* 370 dense.xIonDense[nelem][nelem-ipISO-1]; 371 372 /* now add recombination to next lower stage */ 373 sink += ionbal.RateRecomTot[nelem][nelem-ipISO-1]; 352 374 } 353 375 -
branches/newmole/source/newton_step.cpp
r2069 r2086 93 93 b1vec[i] = ervals[i]; 94 94 } 95 //fprintf(ioQQQ,"Chem "); 95 96 merror = solve_system(amat,b1vec,n); 96 97 error0 = error1; … … 287 288 // fprintf(stdout,"Err %d %g\n",k,maxerr); 288 289 getrs_wrapper('N',n,1,&lufac[0],n,&ipiv[0],&err[0],n,&merror); 290 if (0) 291 { 292 // Quick-and-dirty condition number estimate 293 // see Golub & Van Loan, 3rd Edn, p128 294 double maxb=0., maxe=0.; 295 for (i=0;i<n;i++) 296 { 297 if (fabs(b[i])>maxb) 298 maxb = fabs(b[i]); 299 if (fabs(err[i])>maxe) 300 maxe = fabs(err[i]); 301 } 302 if (k == 0) 303 { 304 fprintf(ioQQQ,"Condition %g\n",1e16*maxe/maxb); 305 } 306 } 289 307 for (i=0;i<n;i++) 290 308 { -
branches/newmole/source/parse_caseb.cpp
r2079 r2086 4 4 #include "cddefines.h" 5 5 #include "opacity.h" 6 #include "rfield.h"6 //#include "rfield.h" 7 7 #include "parse.h" 8 8 … … 19 19 { 20 20 opac.lgCaseB = false; 21 rfield.lgInducProcess = false;21 // rfield.lgInducProcess = false; 22 22 } 23 23 else if( nMatch(" B ",chCard) ) 24 24 { 25 25 opac.lgCaseB = true; 26 rfield.lgInducProcess = false;26 // rfield.lgInducProcess = false; 27 27 } 28 28 else if( nMatch(" C ",chCard) ) -
branches/newmole/source/radius_increment.cpp
r2080 r2086 105 105 * change is to only do this branch if abundance is large enough to be detected by 106 106 * the ionization convergence solvers */ 107 # define ERR_CHK 1.001 108 double err_chk = ERR_CHK; 107 const realnum err_tol = 1e-3; 109 108 /* >>chng 05 sep 02, when low-T solver used solution is approximate, 110 109 * and it must not matter (lot-T solver should not be used if it 111 110 * does matter - so use more liberal expectation */ 112 if( iso.chTypeAtomUsed[ipISO][nelem][0]=='L' ) 113 err_chk *= 10.; 111 //if( iso.chTypeAtomUsed[ipISO][nelem][0]=='L' ) 112 // err_chk *= 10.; 113 realnum abund = 0.; 114 for (long n=0; n<iso.numLevels_local[ipISO][nelem]; ++n) 115 { 116 abund += StatesElem[ipISO][nelem][n].Pop; 117 } 114 118 /* make sure that populations are valid, first check Pop2Ion */ 115 if( StatesElem[ipISO][nelem][0].Pop > 116 dense.xIonDense[nelem][nelem-ipISO]/(double)SDIV(dense.xIonDense[nelem][nelem+1-ipISO])*err_chk && 119 if( iso.pop_ion_ov_neut[ipISO][nelem] > 1e-15 && 120 ! fp_equal_tol(abund*dense.xIonDense[nelem][nelem+1-ipISO],dense.xIonDense[nelem][nelem-ipISO], 121 err_tol*dense.xIonDense[nelem][nelem-ipISO]) 117 122 /*>>chng 06 jul 22, only do check if pops are within tolerance of double */ 118 iso.pop_ion_ov_neut[ipISO][nelem] > 1e-15)123 ) 119 124 { 120 125 fprintf(ioQQQ, … … 122 127 fnzone, 123 128 ipISO , nelem , 124 StatesElem[ipISO][nelem][0].Pop*dense.xIonDense[nelem][nelem+1-ipISO],129 abund*dense.xIonDense[nelem][nelem+1-ipISO], 125 130 dense.xIonDense[nelem][nelem-ipISO] , 126 131 iso.chTypeAtomUsed[ipISO][nelem]); … … 131 136 dense.xIonDense[nelem][nelem+1-ipISO]/SDIV( dense.xIonDense[nelem][nelem-ipISO] ) ); 132 137 fprintf(ioQQQ, 133 "simple %.3e Test shows Pop2Ion (%.6e) >xIonDense[nelem]/[nelem+1] (%.6e) \n",138 "simple %.3e Test shows Pop2Ion (%.6e) != xIonDense[nelem]/[nelem+1] (%.6e) \n", 134 139 iso.xIonSimple[ipISO][nelem], 135 StatesElem[ipISO][nelem][0].Pop,140 abund , 136 141 dense.xIonDense[nelem][nelem-ipISO]/(double)SDIV(dense.xIonDense[nelem][nelem+1-ipISO]) ); 137 142 } … … 142 147 * if ok, rm this and all ref to lgSearch, perhaps conv.h header */ 143 148 ASSERT( !conv.lgSearch ); 144 ASSERT( StatesElem[ipISO][nelem][0].Pop <= 145 dense.xIonDense[nelem][nelem-ipISO]/ 146 (double)SDIV(dense.xIonDense[nelem][nelem+1-ipISO])*err_chk || 147 /*>>chng 06 jul 22, only do check if pops are within tolerance of double */ 148 iso.pop_ion_ov_neut[ipISO][nelem] < 1e-15); 149 150 /* make sure that populations are valid, first check Pop2Ion */ 151 if( StatesElem[ipISO][nelem][0].Pop* 152 dense.xIonDense[nelem][nelem+1-ipISO] > 153 dense.xIonDense[nelem][nelem-ipISO]*err_chk && 154 /*>>chng 06 jul 22, only do check if pops are within tolerance of double */ 155 iso.pop_ion_ov_neut[ipISO][nelem] > 1e-15) 156 { 157 fprintf(ioQQQ," DISASTER PopLo fnz %.2f ipISO %li nelem %li pop_ion_ov_neut %.2e 2pops %.6e %.6e solver:%s\n", 158 fnzone, 159 ipISO , nelem , 160 iso.pop_ion_ov_neut[ipISO][nelem] , 161 StatesElem[ipISO][nelem][0].Pop*dense.xIonDense[nelem][nelem+1-ipISO], 162 dense.xIonDense[nelem][nelem-ipISO]*err_chk , 163 iso.chTypeAtomUsed[ipISO][nelem]); 164 } 165 ASSERT( 166 /*(conv.lgSearch || !conv.nPres2Ioniz) || */ 167 (StatesElem[ipISO][nelem][0].Pop*dense.xIonDense[nelem][nelem+1-ipISO]<= 168 dense.xIonDense[nelem][nelem-ipISO]*err_chk) || 169 /*>>chng 06 jul 22, only do check if pops are within tolerance of double */ 170 iso.pop_ion_ov_neut[ipISO][nelem] < 1e-15); 171 # undef ERR_CHK 149 ASSERT( /*>>chng 06 jul 22, only do check if pops are within tolerance of double */ 150 iso.pop_ion_ov_neut[ipISO][nelem] < 1e-15 || 151 fp_equal_tol(abund*dense.xIonDense[nelem][nelem+1-ipISO],dense.xIonDense[nelem][nelem-ipISO], 152 err_tol*dense.xIonDense[nelem][nelem-ipISO]) ); 172 153 } 173 154 }
