00001
00002
00003
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
00024
00025
00026
00027
00028
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
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
00058 Aul = HydroEinstA( N_(ipLo), N_(ipHi) )*z4;
00059 iso_put_error(ipISO,nelem,ipHi,ipLo,IPRAD,0.001f,0.001f);
00060
00061 ASSERT( Aul > 0.);
00062 }
00063 else
00064 {
00065
00066
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
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
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,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,0.001f);
00098 }
00099 else if( ipLo == 0 && ipHi == 1 )
00100 {
00101
00102 Aul = 8.226*pow((double)(nelem+1.),6.);
00103 iso_put_error(ipISO,nelem,ipHi,ipLo,IPRAD,0.001f,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,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,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,0.001f);
00122 }
00123 }
00124
00125 return (realnum)Aul;
00126 }