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 /*iso_continuum_lower - limit max prin. quan. no. due to continuum lowering processes */ 00004 #include "cddefines.h" 00005 #include "phycon.h" 00006 #include "dense.h" 00007 #include "iso.h" 00008 #include "hydrogenic.h" 00009 #include "trace.h" 00010 00011 void iso_continuum_lower( long ipISO, long nelem ) 00012 { 00013 double a; 00014 long int np, nd, ns, nc; 00015 long eff_charge; 00016 00017 /* size of rate matrices will be defined according to the n calculated here */ 00018 00019 ASSERT( nelem < LIMELM ); 00020 /* this may change at a future date. */ 00021 ASSERT( ipISO <= 1 ); 00022 00023 eff_charge = nelem + 1 - ipISO; 00024 00025 /* Particle packing - the density here should be density of all nuclei in the plasma */ 00026 /* This one is just nuclear charge, which is independent of iso, and always nelem+1. */ 00027 a = sqrt( 1.8887E8 * (nelem+1.) / pow((double)dense.xNucleiTotal, 0.333) ); 00028 ASSERT( a > 0. ); 00029 ASSERT( a < 1.e15 ); 00030 if( a > (double)iso.n_HighestResolved_max[ipISO][nelem]+(double)iso.nCollapsed_max[ipISO][nelem] ) 00031 { 00032 np = iso.n_HighestResolved_max[ipISO][nelem]+iso.nCollapsed_max[ipISO][nelem] + 1; 00033 } 00034 else 00035 np = (long)a; 00036 00037 /* Debye shielding - the density here is electron density */ 00038 /* This one depends on effective charge. */ 00039 a = 2.6E7 * eff_charge * eff_charge * pow( phycon.te/dense.eden, 0.25); 00040 ASSERT( a > 0. ); 00041 ASSERT( a < 1.e15 ); 00042 if( a > (double)iso.n_HighestResolved_max[ipISO][nelem]+(double)iso.nCollapsed_max[ipISO][nelem] ) 00043 { 00044 nd = iso.n_HighestResolved_max[ipISO][nelem]+iso.nCollapsed_max[ipISO][nelem] + 1; 00045 } 00046 else 00047 nd = (long)a; 00048 00049 /* Stark broadening - this should be the density of singly charged ions, 00050 * both positive and negative. The sum of protons, electrons, and HeII should be 00051 * good enough. */ 00052 /* This one depends on effective charge. */ 00053 a = 3171. * pow( (double)eff_charge, 0.8 ) * pow( dense.eden + (double)dense.xIonDense[ipHYDROGEN][1] 00054 + (double)dense.xIonDense[ipHELIUM][1], -0.1333); 00055 ASSERT( a > 0. ); 00056 ASSERT( a < 1.e15 ); 00057 if( a > (double)iso.n_HighestResolved_max[ipISO][nelem]+(double)iso.nCollapsed_max[ipISO][nelem] ) 00058 { 00059 ns = iso.n_HighestResolved_max[ipISO][nelem]+iso.nCollapsed_max[ipISO][nelem] + 1; 00060 } 00061 else 00062 ns = (long)a; 00063 00064 ASSERT( np > 3 ); 00065 ASSERT( nd > 3 ); 00066 ASSERT( ns > 3 ); 00067 00068 nc = MIN3(np, nd, ns); 00069 00070 /* I assert greater than three because the code depends upon having at least up to n=3, and 00071 * because it would take EXTREMELY dense conditions to lower the continuum that much, and 00072 * the code would probably get a wrong answer then anyway. */ 00073 ASSERT( nc > 3 ); 00074 00075 if( nc < iso.n_HighestResolved_max[ipISO][nelem]) 00076 { 00077 iso.lgLevelsLowered[ipISO][nelem] = true; 00078 iso.lgLevelsEverLowered[ipISO][nelem] = true; 00079 hydro.lgReevalRecom = true; 00080 iso.n_HighestResolved_local[ipISO][nelem] = nc; 00081 iso.nCollapsed_local[ipISO][nelem] = 0; 00082 iso.numLevels_local[ipISO][nelem] = iso_get_total_num_levels( ipISO, nc, 0 ); 00083 } 00084 /* Here is the case where the critical n lies among the one or more collapsed levels */ 00085 /* we just get rid of any that are too high. */ 00086 else if( nc <= iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem] ) 00087 { 00088 iso.lgLevelsLowered[ipISO][nelem] = true; 00089 iso.lgLevelsEverLowered[ipISO][nelem] = true; 00090 hydro.lgReevalRecom = true; 00091 iso.n_HighestResolved_local[ipISO][nelem] = iso.n_HighestResolved_max[ipISO][nelem]; 00092 iso.nCollapsed_local[ipISO][nelem] = nc - iso.n_HighestResolved_local[ipISO][nelem]; 00093 iso.numLevels_local[ipISO][nelem] = 00094 iso_get_total_num_levels( ipISO, iso.n_HighestResolved_max[ipISO][nelem], iso.nCollapsed_local[ipISO][nelem] ); 00095 } 00096 /* This is usually where control will flow, because in most conditions the continuum will not be lowered. 00097 * Nothing changes in this case. */ 00098 else 00099 { 00100 iso.numLevels_local[ipISO][nelem] = iso.numLevels_max[ipISO][nelem]; 00101 iso.nCollapsed_local[ipISO][nelem] = iso.nCollapsed_max[ipISO][nelem]; 00102 iso.n_HighestResolved_local[ipISO][nelem] = iso.n_HighestResolved_max[ipISO][nelem]; 00103 00104 /* if levels were lowered on last pass but are not now, must reeval */ 00105 if( iso.lgLevelsLowered[ipISO][nelem] ) 00106 { 00107 hydro.lgReevalRecom = true; 00108 } 00109 else 00110 { 00111 hydro.lgReevalRecom = false; 00112 } 00113 00114 iso.lgLevelsLowered[ipISO][nelem] = false; 00115 } 00116 00117 /* None of these can be greater than that which was originally malloc'd. */ 00118 ASSERT( iso.numLevels_local[ipISO][nelem] <= iso.numLevels_max[ipISO][nelem] ); 00119 ASSERT( iso.nCollapsed_local[ipISO][nelem] <= iso.nCollapsed_max[ipISO][nelem] ); 00120 ASSERT( iso.n_HighestResolved_local[ipISO][nelem] <= iso.n_HighestResolved_max[ipISO][nelem] ); 00121 00122 /* don't print levels in continuum. */ 00123 iso.numPrintLevels[ipISO][nelem] = MIN2( iso.numPrintLevels[ipISO][nelem], iso.numLevels_local[ipISO][nelem] ); 00124 00125 /* Lyman lines can not be greater than original malloc or critical pqn. */ 00126 iso.nLyman[ipISO] = MIN2( nc, iso.nLyman[ipISO]); 00127 00128 if( trace.lgTrace ) 00129 { 00130 fprintf( ioQQQ," iso_continuum_lower: ipISO %li nelem %li nc %li numLevels %li nCollapsed %li n_HighestResolved %li \n", 00131 ipISO, 00132 nelem, 00133 nc, 00134 iso.numLevels_local[ipISO][nelem], 00135 iso.nCollapsed_local[ipISO][nelem], 00136 iso.n_HighestResolved_local[ipISO][nelem] 00137 ); 00138 } 00139 00140 return; 00141 }