firelib.c

这是一个GPS相关的程序

        doEffectiveWind = 1;
        /* But since BEHAVE doesn't calculate effective wind when no spread. */
        /* we wont either. */
        effectiveWind = 0.;
        doEffectiveWind = 0;
        checkWindLimit = 0;
    }

    /* Situation 2: no wind and no wind */
    else if ( phiEw < Smidgen )
    {
        phiEw = 0.;
        effectiveWind = 0.;
        doEffectiveWind = 0;
        spreadMax = Fuel_Spread0(catalog,model);
        azimuthMax = 0;
        checkWindLimit = 0;
    }

    /* Situation 3: wind with no slope. */
    else if ( slope < Smidgen )
    {
        effectiveWind = windFpm;
        doEffectiveWind = 0;
        spreadMax = Fuel_Spread0(catalog,model) * (1. + phiEw);
        azimuthMax = windDeg;
        checkWindLimit = 1;
    }

    /* Situation 4: slope with no wind. */
    else if ( windFpm < Smidgen )
    {
        doEffectiveWind = 1;
        spreadMax = Fuel_Spread0(catalog,model) * (1. + phiEw);
        azimuthMax = upslope;
        checkWindLimit = 1;
    }

    /* Situation 5: wind blows upslope. */
    else if ( Equal(upslope,windDeg) )
    {
        doEffectiveWind = 1;
        spreadMax = Fuel_Spread0(catalog,model) * (1. + phiEw);
        azimuthMax = upslope;
        checkWindLimit = 1;
    }

    /* Situation 6: wind blows cross slope. */
    else
    {
        /* Recalculate spread rate in the optimal direction. */
        splitDeg = (upslope<=windDeg) ? windDeg-upslope : 360.-upslope+windDeg;
        splitRad = DegreesToRadians(splitDeg);
        slpRate  = Fuel_Spread0(catalog,model) * Fuel_PhiSlope(catalog,model);
        wndRate  = Fuel_Spread0(catalog,model) * Fuel_PhiWind(catalog,model);
        x        = slpRate + wndRate * cos(splitRad);
        y        = wndRate * sin(splitRad);
        rv       = sqrt(x*x + y*y);
        spreadMax= Fuel_Spread0(catalog,model) + rv;

        /* Recalculate phiEw in the optimal direction. */
        phiEw    = spreadMax / Fuel_Spread0(catalog,model) - 1.0;
        doEffectiveWind = (phiEw > Smidgen) ? 1 : 0;
        checkWindLimit = 1;

        /* Recalculate direction of maximum spread in azimuth degrees. */
        al = asin(fabs(y) / rv);
        if ( x >= 0. )
            a = (y >= 0.) ? al          : M_PI + M_PI - al;
        else
            a = (y >= 0.) ? (M_PI - al) : (M_PI + al);

        splitDeg = RadiansToDegrees(a);
        if ( (azimuthMax = upslope + splitDeg) > 360. )
            azimuthMax -= 360.;
    }

    /* Recalculate effective wind speed based upon phiEw. */
    if ( doEffectiveWind )
        effectiveWind = pow( (phiEw * Fuel_WindE(catalog,model)),
                             (1. / Fuel_WindB(catalog,model)) );

    /* If effective wind exceeds maximum wind, scale back spread & phiEw. */
    if ( checkWindLimit )
    {
        maxWind = 0.9 * Fuel_RxIntensity(catalog,model);
        if ( effectiveWind > maxWind )
        {
            phiEw = (maxWind < Smidgen) ? 0. :
                Fuel_WindK(catalog,model) * pow(maxWind, Fuel_WindB(catalog,model));

            spreadMax = Fuel_Spread0(catalog,model) * (1. + phiEw);
            effectiveWind = maxWind;
            windLimit = 1;
        }
    }

    /* Determine fire ellipse parameters from the effective wind speed. */
    /* = 1. + 0.25 * (Fuel_EffectiveWind(catalog,model) / 88.0); */
    if ( effectiveWind > Smidgen )
    {
        lwRatio = 1. + 0.002840909 * effectiveWind;
        eccentricity = sqrt(lwRatio * lwRatio - 1.0) / lwRatio;
    }

    /* Store the results. */
    Fuel_Aspect(catalog,model)       = aspect;
    Fuel_WindDir(catalog,model)      = windDeg;
    Fuel_PhiEffWind(catalog,model)   = phiEw;
    Fuel_EffectiveWind(catalog,model)= effectiveWind;
    Fuel_WindLimit(catalog,model)    = windLimit;
    Fuel_SpreadMax(catalog,model)    = Fuel_SpreadAny(catalog,model)  = spreadMax;
    Fuel_AzimuthMax(catalog,model)   = Fuel_AzimuthAny(catalog,model) = azimuthMax;
    Fuel_LwRatio(catalog,model)      = lwRatio;
    Fuel_Eccentricity(catalog,model) = eccentricity;

    /* Initialize behavior variables updated by Fire_SpreadAtAzimuth(). */
    Fuel_ByramsIntensity(catalog,model) = 0.;
    Fuel_FlameLength(catalog,model)     = 0.;
    Fuel_ScorchHeight(catalog,model)    = 0.;

    return (FuelCat_Status(catalog) = FIRE_STATUS_OK);
}

/*
 *******************************************************************************
 *
 *  Fire_SpreadAtAzimuth()
 *
 *  Description
 *      Calculates fire spread rate in a specified direction and optionally
 *      calculates the fireline intensity, flame length, and scorch height
 *      along the fire spread vector.
 *
 *  Side Effects
 *      The following variables are updated:
 *          Fuel_SpreadAny(catalog,model)
 *          Fuel_AzimuthAny(catalog,model) == azimuth;
 *          Fuel_ByramsIntensity(catalog,model) is updated if FIRE_BYRAMS.
 *          Fuel_FlameLength(catalog,model) is updated if FIRE_FLAME.
 *          Fuel_ScorchHeight(catalog,model) is updated if FIRE_SCORCH.
 *  Notes
 *      The calculations depend upon the most recent calls to
 *      Fire_SpreadNoWindNoSlope() and Fire_SpreadWindSlopeMax() for this model.
 *
 *      The input azimuth is the degrees clockwise from north.
 *
 *  Function Returns
 *      FIRE_STATUS_OK or FIRE_STATUS_ERROR.
 *      Return status and error text are stored in the Fire Catalog's buffers.
 *
 *******************************************************************************
 */

int
Fire_SpreadAtAzimuth ( FuelCatalogPtr catalog, size_t model, double azimuth, size_t which )
   // FuelCatalogPtr catalog;     /* FuelCatalogData instance pointer           */
   // size_t  model;              /* fuel model number            [0-maxModels] */
   // double  azimuth;            /* fire spread azimuth     (deg from upslope) */
   // size_t  which;      /* FIRE_NONE | FIRE_BYRAMS | FIRE_FLAME | FIRE_SCORCH */
{
    double dir;
    double radians;
    double byrams;
    double mph;
    size_t lo, hi, mid, n;

    /* Validate catalog and the fuel model. */
    assert(catalog!= NULL && FuelCat_MagicCookie(catalog)==FIRE_CATALOG_MAGIC);
    if ( ! Fire_FuelModelExists(catalog,model) )
    {
        sprintf(FuelCat_Error(catalog),
            "Fire_SpreadAtAzimuth(): el modelo de combustible %d no existe en el cat醠ogo de combustibles \"%s\".",
            model, FuelCat_Name(catalog));
        return (FuelCat_Status(catalog) = FIRE_STATUS_ERROR);
    }

    /* Situation 1: no fire or reaction intensity, so no Byrams or flame. */
    if ( Fuel_SpreadMax(catalog,model) < Smidgen )
        return (FuelCat_Status(catalog) = FIRE_STATUS_OK);

    /* Situation 2: phiEw is zero OR azimuth is in the max spread direction */
    if ( Fuel_PhiEffWind(catalog,model) < Smidgen
      || Equal(Fuel_AzimuthMax(catalog,model),azimuth) )
    {
        Fuel_SpreadAny(catalog,model) = Fuel_SpreadMax(catalog,model);
    }

    /* Situation 3: wind and/or slope and azimuth not in max spread direction */
    else
    {
        /* Angle between maximum spread azimuth and requested azimuth. */
        if ( (dir = fabs(Fuel_AzimuthMax(catalog,model) - azimuth)) > 180. )
            dir = 360. - dir;
        radians = DegreesToRadians(dir);

        /* Calculate the fire spread rate in this azimuth. */
        Fuel_SpreadAny(catalog,model)
            = Fuel_SpreadMax(catalog,model)
            * (1. - Fuel_Eccentricity(catalog,model))
            / (1. - Fuel_Eccentricity(catalog,model) * cos(radians));
    }
    Fuel_AzimuthAny(catalog,model) = azimuth;

    /* Additional fire behavior outputs. */
    if ( which )
    {
        /* Must compute Byram's if any of the three are requested. */
        byrams = Fuel_ResidenceTime(catalog,model)
               * Fuel_SpreadAny(catalog,model)
               * Fuel_RxIntensity(catalog,model)
               / 60.;

        /* Byrams intensity is requested. */
        if ( which & FIRE_BYRAMS )
            Fuel_ByramsIntensity(catalog,model) = byrams;

        /* Flame length is requested. */
        if ( (which & FIRE_FLAME) )
        {
            if ( byrams < Smidgen )
            {
                Fuel_FlameLength(catalog,model) = 0.;
            }
            else
            {
                /* Use lookup table if it exists & includes this intensity. */
                if ( (n = FuelCat_FlameClasses(catalog)) > 0
                  && FuelCat_FlameArray(catalog)[n-1] > byrams )
                {
                    hi = n-1;
                    lo = 0;
                    do {
                        mid = lo + (hi-lo)/2;
                        if ( FuelCat_FlameArray(catalog)[mid] > byrams )
                            hi = mid;
                        else
                            lo = mid + 1;
                    } while (lo != hi);
                    Fuel_FlameLength(catalog,model) =
                        FuelCat_FlameStep(catalog) * (lo+1);
                }
                /* otherwise compute flame length from scratch. */
                else
                {
                    Fuel_FlameLength(catalog,model) = 0.45 * pow(byrams, 0.46);
                }
            }
        }

        /* Scorch height is requested. */
        if ( (which & FIRE_SCORCH) )
        {
            if ( byrams < Smidgen )
            {
                Fuel_ScorchHeight(catalog,model) = 0.;
            }
            else
            {
                mph = Fuel_WindSpeed(catalog,model) / 88.;
                Fuel_ScorchHeight(catalog,model) =
                    pow(byrams, 1.166667) / sqrt(byrams + (mph * mph * mph));
            /*  Fuel_ScorchHeight(catalog,model) *= (63. / (140. - temp_f) ); */
            }
        }
    }

    return (FuelCat_Status(catalog) = FIRE_STATUS_OK);
}

/*
 *******************************************************************************
 *
 *  Fire_FlameScorch()
 *
 *  Description
 *      Calculates the flame length and/or scorch height for the current
 *      Byram's intensity and azimuth (as determined by the most recent
 *      call to Fire_SpreadAtAzimuth()).
 *      Uses the Flame Length Table if it exists.
 *      Offers a method of getting flame length if Fire_SpreadAtAzimuth()
 *      is not calculating it.
 *
 *  Side Effects
 *      The following variables are updated:
 *          Fuel_FlameLength(catalog,model) is updated.
 *          Fuel_ScorchHeight(catalog,model)
 *  Notes
 *      The calculations depend upon the most recent calls to
 *      Fire_SpreadNoWindNoSlope(), Fire_SpreadWindSlopeMax(), and
 *      Fire_SpreadAtAzimuth() for this model.
 *
 *  Function Returns
 *      FIRE_STATUS_OK or FIRE_STATUS_ERROR.
 *      Return status and error text are stored in the Fire Catalog's buffers.
 *
 *******************************************************************************
 */

int
Fire_FlameScorch ( FuelCatalogPtr catalog, size_t model, size_t which )
   // FuelCatalogPtr catalog;     /* FuelCatalogData instance pointer           */
   // size_t  model;              /* fuel model number            [0-maxModels] */
   // size_t  which;      /* FIRE_NONE | FIRE_BYRAMS | FIRE_FLAME | FIRE_SCORCH */
{
    double byrams;
    double mph;
    size_t lo, hi, mid, n;

    /* Validate catalog and the fuel model. */
    assert(catalog!= NULL && FuelCat_MagicCookie(catalog)==FIRE_CATALOG_MAGIC);
    if ( ! Fire_FuelModelExists(catalog,model) )
    {
        sprintf(FuelCat_Error(catalog),
            "Fire_FlameScorch(): el modelo de combustible %d no existe en el cat醠ogo de combustibles \"%s\".",
            model, FuelCat_Name(catalog));
        return (FuelCat_Status(catalog) = FIRE_STATUS_ERROR);
    }

    byrams = Fuel_ResidenceTime(catalog,model)
           * Fuel_SpreadAny(catalog,model)
           * Fuel_RxIntensity(catalog,model)
           / 60.;

    /* Flame length is requested. */
    if ( (which & FIRE_FLAME) )
    {
        if ( byrams < Smidgen )
        {
            Fuel_FlameLength(catalog,model) = 0.;
        }
        else
        {
            /* Use lookup table if it exists & includes this intensity. */
            if ( (n = FuelCat_FlameClasses(catalog)) > 0
              && FuelCat_FlameArray(catalog)[n-1] > byrams )
            {
                hi = n-1;
                lo = 0;
                do {
                    mid = lo + (hi-lo)/2;