/home66/gary/public_html/cloudy/c08_branch/source/dense_fabden.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*dense_fabden called by dlaw command, returns density for any density law */
#include "cddefines.h"
#include "rfield.h"
#include "dense.h"

/*dense_fabden called by dlaw command, returns density for any density law */
double dense_fabden(double radius, 
  double depth)
{
        double a , b;
        a = depth;
        b = radius;

        fprintf(ioQQQ,"PROBLEM DISASTER The default version of dense_fabden is "
                "still in dense_fabden.cpp - edit and replace this file to use your version.\n");
        /* return val is example - must return density */
        return(a*b);

        /* these are two examples of using this routine */
#       if 0
        double fabden_v;
        double z , z0 , density , rad_midplane;
        DEBUG_ENTRY( "dense_fabden()" );

        /* this routine is called when the DLAW command is given,
         * and must return the hydrogen density at the current radius or depth
         * RADIUS is the radius, the distance from the center of symmetry.
         * DEPTH is the depth into the cloud, from the illuminated face
         * both are in cm
         *
         * radius, depth, and the array DensityLaw are double precision, although
         * dense_fabden itself is not
         *
         * this is one way to generate a density
         * dense_fabden = radius * depth
         *
         * this is how to use the parameters on the dlaw command
         * dense_fabden = DensityLaw(1) + radius * DensityLaw(2)
         *
         * following must be removed if this sub is to be used */
        /*fprintf( ioQQQ, " the old version of dense_fabden must be deleted, and a new one put in place.  Sorry.\n" );*/

        if( depth > radius )
                TotalInsanity();

        rad_midplane = radius * cos( dense.DensityLaw[0] );
        z = radius * sin( dense.DensityLaw[0] );
        z0 = 0.047 * pow( rad_midplane/AU , 1.25 ) * AU;
        /* density in gm / cm3 */
        density = 1.4e-9 * pow( rad_midplane/AU , -11./4. ) * sexp( POW2(z/z0) );
        /* convert density to cm-3 */
        fabden_v = density/(ATOMIC_MASS_UNIT * 1.2);

        fprintf(ioQQQ, "bug dense_fabden zone %.2f radius %e density %e \n",
                fnzone , radius , fabden_v );
        return( fabden_v );
#       endif

#if 0

        /* update stromgren radius to most recent value */
        if( rfield.lgUSphON )
                dense.DensityLaw[2] = rfield.rstrom/PARSEC;

        /* this is the density law used in the Wen & O'Dell, 1995, ApJ 438, 784 paper */
        powexp = MIN2((radius/parsec-dense.DensityLaw[2])/dense.DensityLaw[1],dense.DensityLaw[3]);
        fabden_v = pow(10.,dense.DensityLaw[0])*exp(powexp);

        fabden_v = dense.DensityLaw[0];
        /* just here to stop compilers from flagging unused vars */
        temp = radius + depth; 
        if( fabden_v == 0. )
        {
                cdEXIT(EXIT_FAILURE);
        }
        else
        {
                /* when this routine is used the following branch is the correct exit */
                /* following should return correct value of the density at this position, 
                 * temp is there to trick lint */
                return( fabden_v*temp );
        }
#endif
}

---