00001
00002
00003
00004 #include "cddefines.h"
00005 #include "dense.h"
00006 #include "phycon.h"
00007 #include "path.h"
00008 #include "hydrogenic.h"
00009
00010 t_hydrobranch::t_hydrobranch()
00011 {
00012 const long VERSION_MAGIC = 20070109L;
00013 const static char chFile[] = "hydrobranch.dat";
00014 char chString[FILENAME_PATH_LENGTH_2];
00015
00016 DEBUG_ENTRY( "t_hydrobranch()" );
00017
00018 strcpy( chString, ( lgDataPathSet ? chDataPath : "" ) );
00019 strcat( chString, chFile );
00020
00021 FILE *io;
00022 if( (io = fopen(chString,"r")) == NULL )
00023 {
00024 fprintf( ioQQQ, " Could not open %s for reading\n", chString );
00025 fprintf( ioQQQ, " Is the path set correctly?\n" );
00026 puts( "[Stop in t_hydrobranch]" );
00027 cdEXIT(EXIT_FAILURE);
00028 }
00029
00030 long i=-1;
00031
00032 bool lgErr = ( fscanf( io, "%ld", &i ) != 1 );
00033 if( lgErr || i != VERSION_MAGIC )
00034 {
00035 fprintf( ioQQQ, " File %s has incorrect version: %ld\n", chString, i );
00036 fprintf( ioQQQ, " I expected to find version: %ld\n", VERSION_MAGIC );
00037 puts( "[Stop in t_hydrobranch]" );
00038 cdEXIT(EXIT_FAILURE);
00039 }
00040
00041 ReadBlock(io,Alow);
00042 ReadBlock(io,Blow);
00043 ReadBlock(io,Clow);
00044 ReadBlock(io,Ahigh);
00045 ReadBlock(io,Aint);
00046 ReadBlock(io,Bint);
00047 ReadBlock(io,Cint);
00048 ReadBlock(io,Azero);
00049 ReadBlock(io,Bzero);
00050 ReadBlock(io,Czero);
00051
00052
00053
00054 int nd2=-1;
00055
00056 lgErr = lgErr || ( fscanf( io, "%i", &nd2 ) != 1 );
00057 ASSERT( !lgErr && nd2 == NDIM2 );
00058
00059 for( int j=0; j < nd2; j++ )
00060 lgErr = lgErr || ( fscanf( io, "%le", &DenCut[j] ) != 1 );
00061
00062 ASSERT( !lgErr );
00063
00064 DEBUG_EXIT( "t_hydrobranch()" );
00065 }
00066
00067 void t_hydrobranch::ReadBlock(FILE* io,
00068 double x[][NDIM2])
00069 {
00070 int nd1, nd2;
00071
00072 DEBUG_ENTRY( "ReadBlock()" );
00073
00074 bool lgErr = ( fscanf( io, "%i %i", &nd1, &nd2 ) != 2 );
00075
00076 ASSERT( !lgErr && nd1 == NDIM1 && nd2 == NDIM2 );
00077
00078 for( int i=0; i < nd1; i++ )
00079 for( int j=0; j < nd2; j++ )
00080 lgErr = lgErr || ( fscanf( io, "%le", &x[i][j] ) != 1 );
00081
00082 ASSERT( !lgErr );
00083
00084 DEBUG_EXIT( "ReadBlock()" );
00085 return;
00086 }
00087
00088 double t_hydrobranch::HydroBranchFunc(
00089
00090 long int nHi,
00091
00092 long int nLo,
00093
00094 long int ipZPhys) const
00095 {
00096
00097
00098
00099
00100
00101 long int nelem,
00102 ipHi,
00103 ipLo;
00104 double DenScale,
00105 RedMass,
00106 answer,
00107 den,
00108 denhigh,
00109 denint,
00110 denlow,
00111 denzero,
00112 t,
00113 teuse,
00114 tt,
00115 z7;
00116
00117 DEBUG_ENTRY( "HydroBranchFunc()" );
00118
00119 ASSERT( dense.eden > 0. );
00120
00121
00122 ipLo = MIN2( nLo , nHi ) - 1;
00123 ipHi = MAX2( nLo , nHi ) - 1;
00124
00125
00126 ASSERT(ipLo >= 1 );
00127
00128 ASSERT( ipHi >= 3 );
00129
00130 ASSERT( ipHi <= 14);
00131
00132
00133 ASSERT( ipZPhys > 0 );
00134
00135 nelem = ipZPhys - 1;
00136
00137 RedMass = dense.AtomicWeight[nelem]*1./
00138 (dense.AtomicWeight[nelem] + 1.);
00139
00140 z7 = powi((double)ipZPhys,7);
00141
00142 DenScale = log10(sqrt(2.*RedMass)/z7);
00143
00144 den = log10(dense.eden) + DenScale;
00145
00146
00147 den = MAX2( -2. , den );
00148
00149
00150
00151
00152 t = phycon.te/POW2((double)ipZPhys);
00153 teuse = MAX2(100.,t);
00154 teuse = MIN2(3e4,teuse);
00155 t = log10(teuse);
00156
00157 tt = phycon.te/POW2((double)ipZPhys);
00158 teuse = MAX2(500.,tt);
00159 teuse = MIN2(3e4,teuse);
00160 tt = log10(teuse);
00161
00162
00163 denlow = Alow[ipLo][ipHi] + t*(Blow[ipLo][ipHi] + t*Clow[ipLo][ipHi]);
00164
00165 denhigh = Ahigh[ipLo][ipHi];
00166
00167
00168 denzero = Azero[ipLo][ipHi] + t*(Bzero[ipLo][ipHi] + t*Czero[ipLo][ipHi]);
00169
00170
00171 denint = Aint[ipLo][ipHi] + tt*(Bint[ipLo][ipHi] + tt*Cint[ipLo][ipHi]);
00172
00173
00174 if( den >= (DenCut[ipHi] + 2.) )
00175 {
00176
00177 answer = denhigh;
00178 }
00179 else if( den >= DenCut[ipHi] )
00180 {
00181
00182 answer = denint + (denhigh - denint)*(den - DenCut[ipHi])/2.0;
00183 }
00184 else if( den <= 4.0 )
00185 {
00186
00187 if( ipZPhys != 1 && t <= 3.3 )
00188 {
00189
00190 answer = denzero + (denlow - denzero)*(den - 2.0 - DenScale)/(2.0 - DenScale);
00191 if( answer <= 0. )
00192 {
00193 fprintf(ioQQQ," PROBLEM hydrobranch low den limit ans %g %g %g %g %g\n",
00194 answer , denzero , denlow , den , DenScale);
00195 TotalInsanity();
00196 }
00197 }
00198 else
00199 {
00200 answer = denlow;
00201 }
00202 }
00203 else
00204 {
00205
00206 answer = denlow + (denint - denlow)*(den - 4.0)/5.0;
00207 }
00208
00209 if( answer <= 0. )
00210 {
00211 fprintf(ioQQQ,
00212 " HydroBranchFunc finds insane branching ratio, Z=%ld lo=%ld up=%ld val=%g\n" ,
00213 ipZPhys, ipLo+1, ipHi+1, answer);
00214 puts( "[Stop in HydroBranchFunc]" );
00215 cdEXIT(EXIT_FAILURE);
00216 }
00217
00218 DEBUG_EXIT( "HydroBranchFunc()" );
00219 return answer;
00220 }