/home66/gary/public_html/cloudy/c08_branch/source/parse_radius.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 */
/*ParseRadius parse the radius command */
#include "cddefines.h"
/*#define       PARSCL  18.489396*/
#include "physconst.h"
#include "optimize.h"
#include "radius.h"
#include "iterations.h"
#include "input.h"
#include "parse.h"

void ParseRadius(char *chCard, 
                 realnum *ar1)
{
        bool lgEOL, 
          lgR2set, 
          lgRLog;
        long int i;
        double a, 
          convl;

        DEBUG_ENTRY( "ParseRadius()" );

        /* log of inner and outer radii, default second=infinity,
         * if R2<R1 then R2=R1+R2
         * there is an optional keyword, "PARSEC" on the line, to use PC as units */
        if( nMatch("PARS",chCard) )
        {
                /*>>chng 06 mar 18, from above to below */
                /*convl = PARSCL;*/
                convl = log10( PARSEC );
        }
        else
        {
                convl = 0.;
        }

        /* if linear appears on line, then radius is linear, otherwise, log */
        if( nMatch("LINE",chCard) )
        {
                lgRLog = false;
        }
        else
        {
                lgRLog = true;
        }

        i = 5;
        *ar1 = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
        if( lgEOL )
        {
                fprintf( ioQQQ, " There should have been a number on this line.  Sorry.\n" );
                cdEXIT(EXIT_FAILURE);
        }

        /* option for linear or log radius */
        if( lgRLog )
        {
                *ar1 += (realnum)convl;
        }
        else
        {
                if( *ar1 > 0. )
                {
                        *ar1 = (realnum)(log10(*ar1) + convl);
                }
                else
                {
                        fprintf(ioQQQ,"The first radius is negative and linear radius is set - this is impossible.\n");
                        cdEXIT(EXIT_FAILURE);
                }
        }

        if( *ar1 > 37 || *ar1 < -37. )
        {
                fprintf(ioQQQ,"WARNING - the log of the radius is %e - this is too big and I would crash.\n", *ar1 );
                fprintf(ioQQQ," Sorry.\n" );
                cdEXIT(EXIT_FAILURE);
        }

        radius.Radius = pow((realnum)10.f ,*ar1);
        radius.lgRadiusKnown = true;

        /* check for second number, which indicates thickness or outer radius of model */
        a = (double)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
        if( lgEOL )
        {
                /* not set */
                lgR2set = false;
        }
        else
        {
                /* outer radius is set, */
                lgR2set = true;

                /* log or linear option is still in place */
                if( lgRLog )
                {
                        a += convl;
                }

                else
                {
                        /* linear radius - convert to log but first make sure that a is > 0 */
                        if( a > 0. )
                        {
                                a = log10(a) + convl;
                        }
                        else
                        {
                                fprintf(ioQQQ,"The second radius is negative and linear radius is set - this is impossible.\n");
                                cdEXIT(EXIT_FAILURE);
                        }
                }

                if( a > 37 || a < -37. )
                {
                        fprintf(ioQQQ,"WARNING - the log of the outer radius is %e - this is too big and I will soon crash.\n", a );
                        /* flush buffers since we shall soon throw an fpe */
                        fflush( ioQQQ );
                }
                a = pow(10.,a);
                /* check whether it was thickness or outer radius,
                 * we want router to be total thickness of modeled region,
                 * NOT outer radius */
                if( a > radius.Radius )
                {
                        radius.router[0] = a - radius.Radius;
                }

                else
                {
                        radius.router[0] = a;
                }

                for( i=1; i < iterations.iter_malloc; i++ )
                {
                        radius.router[i] = radius.router[0];
                }
        }

        /* vary option */
        if( optimize.lgVarOn )
        {
                /* pointer to where to write */
                optimize.nvfpnt[optimize.nparm] = input.nRead;
                optimize.vincr[optimize.nparm] = 0.5;

                /* flag saying second outer radius or thickness was set */
                if( lgR2set )
                {
                        strcpy( optimize.chVarFmt[optimize.nparm], "RADIUS %f depth or outer R %f" );
                        optimize.nvarxt[optimize.nparm] = 2;
                        /* second number is thickness or outer radius */
                        optimize.vparm[1][optimize.nparm] = (realnum)log10(radius.router[0]);
                        fprintf(ioQQQ,
                                " WARNING - outer radius or thickness was set with a variable radius.\n");
                        fprintf(ioQQQ,
                                " The interpretation of the second number can change from radius to depth as radius changes.\n");
                        fprintf(ioQQQ,
                                " Do not use the second parameter unless you are quite certain that you know what you are doing.\n");
                        fprintf(ioQQQ,
                                " Consider using the STOP THICKNESS command instead.\n");
                }
                else
                {
                        strcpy( optimize.chVarFmt[optimize.nparm], "RADIUS= %f" );
                        optimize.nvarxt[optimize.nparm] = 1;
                }

                /* log of radius is first number */
                optimize.vparm[0][optimize.nparm] = (realnum)log10(radius.Radius);
                ++optimize.nparm;
        }
        return;
}

---