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 /*HydroEinstA calculates Einstein A's from osillator strengths*/ 00004 #include "cddefines.h" 00005 #include "hydro_bauman.h" 00006 #include "hydrooscilstr.h" 00007 #include "hydroeinsta.h" 00008 #include "iso.h" 00009 #include "physconst.h" 00010 #include "taulines.h" 00011 00012 double HydroEinstA(long int n1, 00013 long int n2) 00014 { 00015 long int lower, iupper; 00016 double EinstA_v, 00017 ryd, 00018 xl, 00019 xmicron, 00020 xu; 00021 00022 DEBUG_ENTRY( "HydroEinstA()" ); 00023 /* (lower,upper) of Johnson 1972. */ 00024 00025 /* strictly n -> n' transition probabilities 00026 * no attempt to distribute according to l,l' */ 00027 00028 /* sort out the order of upper and lower, so can be called either way */ 00029 lower = MIN2( n1 , n2 ); 00030 iupper = MAX2( n1, n2 ); 00031 if( lower < 1 || lower == iupper ) 00032 { 00033 fprintf(ioQQQ," HydroEinstA called with impossible ns, =%li %li\n", lower, iupper); 00034 cdEXIT(EXIT_FAILURE); 00035 } 00036 00037 xl = (double)lower; 00038 xu = (double)iupper; 00039 ryd = 1./POW2(xl) - 1./POW2(xu); 00040 xmicron = 1.E4/(ryd*RYD_INF); 00041 EinstA_v = HydroOscilStr(xl,xu)*TRANS_PROB_CONST*1e8f/(POW2(xmicron))*xl*xl/xu/xu; 00042 return( EinstA_v ); 00043 } 00044 00045 realnum hydro_transprob( long nelem, long ipHi, long ipLo ) 00046 { 00047 double Aul, Aul1; 00048 long ipISO = ipH_LIKE; 00049 /* charge to 4th power, needed for scaling laws for As*/ 00050 double z4 = POW4((double)nelem+1.); 00051 DEBUG_ENTRY( "hydro_transprob()" ); 00052 00053 if( ipHi >= iso.numLevels_max[ipISO][nelem]-iso.nCollapsed_max[ipISO][nelem] ) 00054 { 00055 if( ipLo >= iso.numLevels_max[ipISO][nelem]-iso.nCollapsed_max[ipISO][nelem] ) 00056 { 00057 /* Neither upper nor lower is resolved...Aul's are easy. */ 00058 Aul = HydroEinstA( N_(ipLo), N_(ipHi) )*z4; 00059 iso_put_error(ipISO,nelem,ipHi,ipLo,IPRAD,0.001f); 00060 00061 ASSERT( Aul > 0.); 00062 } 00063 else 00064 { 00065 /* Lower level resolved, upper not. First calculate Aul 00066 * from upper level with ang mom one higher. */ 00067 Aul = H_Einstein_A( N_(ipHi), L_(ipLo)+1, N_(ipLo), L_(ipLo), nelem+1 ); 00068 00069 iso.CachedAs[ipISO][nelem][ N_(ipHi)-iso.n_HighestResolved_max[ipISO][nelem]-1 ][ ipLo ][0] = (realnum)Aul; 00070 00071 Aul *= (2.*L_(ipLo)+3.) * 2. / (2.*(double)N_(ipHi)*(double)N_(ipHi)); 00072 00073 if( L_(ipLo) != 0 ) 00074 { 00075 /* For all l>0, add in transitions from upper level with ang mom one lower. */ 00076 Aul1 = H_Einstein_A( N_(ipHi), L_(ipLo)-1, N_(ipLo), L_(ipLo), nelem+1 ); 00077 00078 iso.CachedAs[ipISO][nelem][ N_(ipHi)-iso.n_HighestResolved_max[ipISO][nelem]-1 ][ ipLo ][1] = (realnum)Aul1; 00079 00080 /* now add in this part after multiplying by stat weight for lHi = lLo-1. */ 00081 Aul += Aul1*(2.*L_(ipLo)-1.) * 2. / (2.*(double)N_(ipHi)*(double)N_(ipHi)); 00082 } 00083 else 00084 iso.CachedAs[ipISO][nelem][ N_(ipHi)-iso.n_HighestResolved_max[ipISO][nelem]-1 ][ ipLo ][1] = 0.f; 00085 00086 iso_put_error(ipISO,nelem,ipHi,ipLo,IPRAD,0.01f); 00087 ASSERT( Aul > 0.); 00088 } 00089 } 00090 else 00091 { 00092 if( N_(ipHi) == N_(ipLo) ) 00093 { 00096 Aul = SMALLFLOAT; 00097 iso_put_error(ipISO,nelem,ipHi,ipLo,IPRAD,0.001f); 00098 } 00099 else if( ipLo == 0 && ipHi == 1 ) 00100 { 00101 /* 2s two photon */ 00102 Aul = 8.226*pow((double)(nelem+1.),6.); 00103 iso_put_error(ipISO,nelem,ipHi,ipLo,IPRAD,0.001f); 00104 } 00105 else if( ipLo == 0 && ipHi == 2 ) 00106 { 00107 Aul = 6.265e8*z4; 00108 iso_put_error(ipISO,nelem,ipHi,ipLo,IPRAD,0.001f); 00109 } 00110 else if( abs( L_(ipLo) - L_(ipHi) )== 1 ) 00111 { 00112 Aul = H_Einstein_A( N_(ipHi), L_(ipHi), N_(ipLo), L_(ipLo), nelem+1 ); 00113 iso_put_error(ipISO,nelem,ipHi,ipLo,IPRAD,0.001f); 00114 } 00115 else 00116 { 00117 ASSERT( N_(ipHi) > N_(ipLo) ); 00118 ASSERT( (L_(ipHi) == L_(ipLo)) || 00119 ( abs(L_(ipHi)-L_(ipLo)) > 1) ); 00120 Aul = SMALLFLOAT; 00121 iso_put_error(ipISO,nelem,ipHi,ipLo,IPRAD,0.001f); 00122 } 00123 } 00124 00125 return (realnum)Aul; 00126 }