00001
00002
00003
00004 #include "cddefines.h"
00005 #include "rfield.h"
00006 #include "dense.h"
00007
00008 double dense_fabden(double radius,
00009 double depth)
00010 {
00011 double a , b;
00012 a = depth;
00013 b = radius;
00014
00015 fprintf(ioQQQ,"The default version of dense_fabden is still in dense_fabden.c - to use this option replace with your version.\n");
00016
00017 return(a*b);
00018
00019 # if 0
00020 double fabden_v;
00021 double z , z0 , density , rad_midplane;
00022 DEBUG_ENTRY( "dense_fabden()" );
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 if( depth > radius )
00043 TotalInsanity();
00044
00045 rad_midplane = radius * cos( dense.DensityLaw[0] );
00046 z = radius * sin( dense.DensityLaw[0] );
00047 z0 = 0.047 * pow( rad_midplane/AU , 1.25 ) * AU;
00048
00049 density = 1.4e-9 * pow( rad_midplane/AU , -11./4. ) * sexp( POW2(z/z0) );
00050
00051 fabden_v = density/(ATOMIC_MASS_UNIT * 1.2);
00052
00053 fprintf(ioQQQ, "bug dense_fabden zone %.2f radius %e density %e \n",
00054 fnzone , radius , fabden_v );
00055
00056 DEBUG_EXIT( "dense_fabden()" );
00057 return( fabden_v );
00058 # endif
00059
00060 #if 0
00061
00062
00063 if( rfield.lgUSphON )
00064 dense.DensityLaw[2] = rfield.rstrom/PARSEC;
00065
00066
00067 powexp = MIN2((radius/parsec-dense.DensityLaw[2])/dense.DensityLaw[1],dense.DensityLaw[3]);
00068 fabden_v = pow(10.,dense.DensityLaw[0])*exp(powexp);
00069
00070 fabden_v = dense.DensityLaw[0];
00071
00072 temp = radius + depth;
00073 if( fabden_v == 0. )
00074 {
00075 puts( "[Stop in dense_fabden]" );
00076 cdEXIT(EXIT_FAILURE);
00077 }
00078 else
00079 {
00080
00081
00082 DEBUG_EXIT( "dense_fabden()" );
00083
00084
00085 return( fabden_v*temp );
00086 }
00087 #endif
00088 }
00089