#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>
#include "fbp97.h"

/* int Date2Julian(); */

FILE *infile, *outfile, *intfile;

/********************** Kerry Anderson's Stuff ***************************/

/* FBP Structure is stored in fbp97.h

   Typical call from C:  

   FBP fbp;   (declaration of the fbp variable structure)
   status = FBPCalc(&fbp);

   Discrepancies between BMW's fbp.c and ST-X-3 

   eqn 32: BMW has 33.5 instead of 35.5 (see ROScalc and Slopecalc)
   eqn 57: BMW had forced SFC to 2.0 for C6 (this is now gone)
   b[O1b]: BMW has 0.0829 instead of 0.0310
   c[O1a]: BMW has 1.41 instead of 1.4
*/

// Note: There is no attempt to account for greened-up deciduous (D2)

int FBPCalc (FBP *fbp)
{
/* inputs */
char   FuelType[4]; 
int    Accel,           /* 1 = point, 0 = line (no longer accepts 2 as a line source) */
       Dj,              /* Julian Day */
       Do,              /* Julian day of minimum FMC */
       ELV,             /* Elevation (m ASL) */
       BUIEff;          /* BUI effect: 1 = yes, 0 = no */
double t,               /* Hours since ignition */
       FFMC,            /* FFMC */
       ISI,             /* ISI */
       BUI,             /* BUI */
       WS,              /* wind speed (kmh) */
       WD,              /* wind direction (degrees) */
       GS,              /* Slope (percent) */
       Aspect,          /* Aspect (degrees) */
       PC,              /* Percent Confier for M1/M2 */
       PDF,             /* Percent Dead Fir for M3/M4 */
       C,               /* Percent Cured for O1a/O1b (85% default) */
       GFL,             /* Grass Fuel Load (kg/m^2) (0.3 default) */
       CBH,             /* Crown to Base Height (FBP defaults)*/
       LAT,             /* Latitude (decimal degrees) */
       LON,             /* Longitude (decimal degrees) */
       FMC,             /* FMC if known */
       theta;           /* elliptical direction of calculation w.r.t. North (degrees)*/

/* outputs */
double WAZ, SAZ, SFC, WSV, RAZ, 
       ROS, FROS, BROS, TROS, /* head, flank, back, theta */
	   ROSt, FROSt, BROSt, TROSt, 
	   CFB, FCFB, BCFB, TCFB, 
	   TFC, FTFC, BTFC, TTFC,
	   HFI, FFI, BFI, TFI,  
	   TI, FTI, BTI, TTI, 
	   LB, E;

int Fuel;

   strcpy(FuelType, fbp->FuelType); 

   Accel  = fbp->Accel;
   Dj     = fbp->Dj;
   Do     = fbp->Do;
   ELV    = fbp->ELV;
   BUIEff = fbp->BUIEff;
   t      = fbp->t;
   FFMC   = fbp->FFMC;
   ISI    = fbp->ISI;
   BUI    = fbp->BUI;
   WS     = fbp->WS;
   WD     = fbp->WD * PI/180.;
   GS     = fbp->GS;
   Aspect = fbp->Aspect * PI/180.;
   PC     = fbp->PC;
   PDF    = fbp->PDF;
   C      = fbp->C;
   GFL    = fbp->GFL;
   CBH    = fbp->CBH;
   LAT    = fbp->LAT;
   LON    = fbp->LON;
   FMC    = fbp->FMC;
   theta  = fbp->theta * PI/180.;
		  
/* Some default values */
      if (Accel < 0) Accel = 0; /* line (no accelleration effect) */
//      if (t < 0.0) t = 0.0;
      if (t < 0.0) t = -t;
      if (FFMC < 0.0) FFMC = 0.0;
      if (BUI < 0.0) BUI = 0.0;
      if (WS < 0.0) WS = 0.0;
      if (GS < 0.0) GS = 0.0;
      if (PC < 0.0) PC = 50.0;
      if (PDF < 0.0) PDF = 35.0;
      if (C <= 0.0) C = 95.;
      if (GFL <= 0.0) GFL = 0.3;
      if (CBH < 0.0) CBH = 7.0;
      if (Dj < 0) Dj = 0;
      if (Do < 0) Do = 0;

/* Convert time from hours to minutes */
   t = t*60.;

/* Corrections to reorient WAZ, SAZ */
   WAZ = WD + PI;
   if (WAZ > 2*PI) WAZ = WAZ - 2*PI;

/* nb: BMW's data set appears to have aspect not saz */  
   SAZ = Aspect + PI;
   if (SAZ > 2*PI) SAZ = SAZ - 2*PI;

/* Make LON positive for the Western Hemisphere */
   LON = -LON;

/* Enter FBP Calculations */
   Fuel = FBPFuel(FuelType);
   if(Fuel==-1 || FFMC == 0.)
   {
      fbp->ROS   = 0.;
      fbp->FROS  = 0.;
      fbp->BROS  = 0.;
      fbp->TROS  = 0.;
      fbp->HROSt = 0.;
      fbp->FROSt = 0.;
      fbp->BROSt = 0.;
      fbp->TROSt = 0.;
      fbp->CFB   = 0.;
      fbp->FCFB  = 0.;
      fbp->BCFB  = 0.;
      fbp->TCFB  = 0.;
      fbp->HFI   = 0.;
      fbp->FFI   = 0.;
      fbp->BFI   = 0.;
      fbp->TFI   = 0.;
      fbp->TFC   = 0.;
      fbp->FTFC  = 0.;
      fbp->BTFC  = 0.;
      fbp->TTFC  = 0.;
      fbp->SFC   = 0.;
      fbp->TI    = -999.;
      fbp->FTI   = -999.;
      fbp->BTI   = -999.;
      fbp->TTI   = -999.;
      fbp->LB    = -999.;
      fbp->RAZ    = -999.;
      fbp->WSV    = -999.;

      return(0);
   }

   CFB=0.;

   if (FMC <= 0)
	  FMC = FMCcalc(LAT, LON, ELV, Dj, Do);

   SFC = SFCcalc(Fuel, FFMC, BUI, PC, GFL);

   if (BUIEff != 1) BUI = 0.0;  /* This turns off BUI effect */

   if (GS > 0. && FFMC > 0.)
   {
      Slopecalc(Fuel, FFMC, BUI, WS, WAZ, GS, SAZ, FMC, SFC, PC, PDF, C, &RAZ, &WSV);
   }
   else
   {
      WSV = WS;
      RAZ = WAZ;
   }

   if (FFMC > 0.)
      ISI = ISICalc(FFMC, WSV);

   ROS = ROScalc(Fuel, ISI, BUI, FMC, SFC, PC, PDF, C);

   LB = LBcalc(Fuel, WSV);
   BROS = BROScalc(Fuel, FFMC, BUI, WSV, FMC, SFC, PC, PDF, C);
   FROS = FROScalc(ROS, BROS, LB);

/* TROS is the rate of spread towards angle theta */
   E = sqrt(1.-1./LB/LB);      /* eccentricity */
   TROS= ROS * (1.-E)/
          (1.-E*cos(theta - RAZ));

   if(Accel == 0)
   {
      ROSt  = ROS;
      FROSt = FROS;
      BROSt = BROS;
      TROSt = TROS;
   }
   else
   {
      ROSt  = ROStcalc(Fuel, ROS, t, CFB);
      FROSt = ROStcalc(Fuel, FROS, t, CFB);
      BROSt = ROStcalc(Fuel, BROS, t, CFB);
      TROSt = ROStcalc(Fuel, TROS, t, CFB);
   }

   if (CFL[Fuel] == 0.)
   {
      CFB = 0.;
      FCFB = 0.;
      BCFB = 0.;
      TCFB = 0.;
   }
   else
   {
      if (!strcasecmp(Fuels[Fuel], "C6"))  
      {
         FCFB = 0.;
         BCFB = 0.;
         TCFB = 0.;
      }
      else
      {
         CFB  = CFBcalc(Fuel, FMC, SFC, ROSt);
         FCFB = CFBcalc(Fuel, FMC, SFC, FROSt);
         BCFB = CFBcalc(Fuel, FMC, SFC, BROSt);
         TCFB = CFBcalc(Fuel, FMC, SFC, TROSt);
      }
   }

   TFC  = TFCcalc(Fuel, CFB, SFC, PC);
   FTFC = TFCcalc(Fuel, FCFB, SFC, PC);
   BTFC = TFCcalc(Fuel, BCFB, SFC, PC);
   TTFC = TFCcalc(Fuel, TCFB, SFC, PC);

   HFI = FIcalc(TFC,  ROSt);
   FFI = FIcalc(FTFC, FROSt);
   BFI = FIcalc(BTFC, BROSt);
   TFI = FIcalc(TTFC, TROSt);

/* For now... */
   TI = 0.;
   FTI = 0.;
   BTI = 0.;
   TTI = 0.;

   fbp->ROS   = ROS;
   fbp->FROS  = FROS;
   fbp->BROS  = BROS;
   fbp->TROS  = TROS;

   if (fbp->t < 0)
   {
		fbp->HROSt = -ROSt;
		fbp->FROSt = -FROSt;
		fbp->BROSt = -BROSt;
		fbp->TROSt = -TROSt;
		fbp->CFB   = -CFB;
   }
   else
   {
		fbp->HROSt = ROSt;
		fbp->FROSt = FROSt;
		fbp->BROSt = BROSt;
		fbp->TROSt = TROSt;
   }

   fbp->FCFB  = FCFB;
   fbp->BCFB  = BCFB;
   fbp->TCFB  = TCFB;
   fbp->HFI   = HFI;
   fbp->FFI   = FFI;
   fbp->BFI   = BFI;
   fbp->TFI   = TFI;
   fbp->TFC   = TFC;
   fbp->FTFC  = FTFC;
   fbp->BTFC  = BTFC;
   fbp->TTFC  = TTFC;
   fbp->SFC   = SFC;
   fbp->TI    = TI;
   fbp->FTI   = FTI;
   fbp->BTI   = BTI;
   fbp->TTI   = TTI;
   fbp->LB    = LB;
   fbp->RAZ   = RAZ * 180./PI;
   fbp->WSV   = WSV;

   return(1);
}

/* Determine and index for the FBP fuel type used throughout the program 
   Note that the order can be adjusted by varying the order in the static variables */

int FBPFuel(char * FuelType)
{
int i, Fuel;
char temp[4];

   if (!strcasecmp(FuelType, "O1"))
      strcpy(temp, "O1b");
   else
	  strcpy(temp, FuelType);

   Fuel = -1;
   for (i = 0; i < MAX_FUELS; i++)
      if (!strcasecmp(FuelType, Fuels[i]) && strcasecmp(FuelType, "WA") && strcasecmp(FuelType, "NF")) Fuel = i;

   return(Fuel);
}

/* Foliar Moisture Content (FMC) calculation 
   Note that 0.5 is added before the integer conversion in equations 2 and 4
   Note that equations 1 and 3 use positive longitude values for Canada */

double FMCcalc(double LAT, double LON, int ELV, int Dj, int Do)
  {
    double FMC, LATN;
    int ND;

    FMC = -1.;

   if (Do <= 0)      /* if Do, date of min FMC, is not known then Do <= 0 */
   {
      if (ELV <= 0)
      {
         LATN = 46.0 + 23.4 * exp(-0.0360 *(150. - LON));               /* 1 */
         Do = (int)(151. * LAT/LATN + 0.5);                      /* 2 (+0.5) */
      }
      else
      {
         LATN = 43. + 33.7 * exp(-0.0351*(150 - LON));                  /* 3 */
         Do = (int)(142.1 * LAT/LATN + (0.0172 * ELV) + 0.5);    /* 4 (+0.5) */
      }
   }

   ND = abs(Dj - Do);                                                   /* 5 */
   
   if (ND < 30)
      FMC = 85. + 0.0189 * ND * ND;                                     /* 6 */
   else 
      if (ND < 50)
         FMC = 32.9 + 3.17 * ND - 0.0288 * ND * ND;                     /* 7 */
      else
         FMC  = 120.;                                                   /* 8 */

   return(FMC);
}

/* Surface Fuel Consumption (SFC) calculation */

double SFCcalc(int Fuel, double FFMC, double BUI, double PC, double GFL)
{
double SFC, FFC, WFC;

   SFC = -1.;

   if (!strcasecmp(Fuels[Fuel], "C1"))
      SFC = 1.5 * (1. - exp(-0.230 * (FFMC - 81)));                    /* 9 */

   if (!strcasecmp(Fuels[Fuel], "C2") ||
       !strcasecmp(Fuels[Fuel], "M3") ||
       !strcasecmp(Fuels[Fuel], "M4"))
      SFC = 5.0 * (1. - exp(-0.0115 * BUI));                          /* 10 */

   if (!strcasecmp(Fuels[Fuel], "C3") ||
       !strcasecmp(Fuels[Fuel], "C4"))
      SFC = 5.0 * pow(1. - exp(-0.0164*BUI), 2.24);                    /* 11 */

   if (!strcasecmp(Fuels[Fuel], "C5") ||
       !strcasecmp(Fuels[Fuel], "C6"))
      SFC = 5.0 * pow(1. - exp(-0.0149*BUI), 2.48);                    /* 12 */

   if (!strcasecmp(Fuels[Fuel], "C7"))
   {
      FFC = 2. * (1. - exp(-0.104 * (FFMC - 70)));                     /* 13 */
      WFC = 1.5 * (1. - exp(-0.0201 * BUI));                           /* 14 */
      SFC = FFC + WFC;                                                 /* 15 */
   }

   if (!strcasecmp(Fuels[Fuel], "D1"))
      SFC = 1.5 * (1. - exp(-0.0183 * BUI));                           /* 16 */

   if (!strcasecmp(Fuels[Fuel], "M1") ||
       !strcasecmp(Fuels[Fuel], "M2"))
   {
      SFC = (PC/100. * SFCcalc(FBPFuel("C2"), FFMC, BUI, PC, GFL))
      + ((100. - PC)/100*SFCcalc(FBPFuel("D1"), FFMC, BUI, PC, GFL));  /* 17 */
   }

   if (!strcasecmp(Fuels[Fuel], "O1a") ||
       !strcasecmp(Fuels[Fuel], "O1b"))
      SFC = GFL;                                                       /* 18 */

   if (!strcasecmp(Fuels[Fuel], "S1"))
   {
      FFC = 4.0 * (1. - exp(-0.025*BUI));                              /* 19 */
      WFC = 4.0 * (1. - exp(-0.034*BUI));                              /* 20 */
      SFC = FFC + WFC;                                                 /* 25 */
   }

   if (!strcasecmp(Fuels[Fuel], "S2"))
   {
      FFC = 10.0 * (1. - exp(-0.013*BUI));                             /* 19 */
      WFC = 6.0 * (1. - exp(-0.060*BUI));                              /* 20 */
      SFC = FFC + WFC;                                                 /* 25 */
   }

   if (!strcasecmp(Fuels[Fuel], "S3"))
   {
      FFC = 12.0 * (1. - exp(-0.0166*BUI));                            /* 19 */
      WFC = 20.0 * (1. - exp(-0.0210*BUI));                            /* 20 */
      SFC = FFC + WFC;                                                 /* 25 */
   }
   if (SFC <= 0.) SFC = 0.000001;
   return(SFC);
}

/* Rate of Spread calculations */

double ROScalc(int Fuel, double ISI, double BUI, double FMC, double SFC, double PC, double PDF, double C)
{
double RSI, CF, BE, ROS;


   RSI = -1.;

/* Note that only preliminary RSS calculations are done for C6 in this routine */
   if (!strcasecmp(Fuels[Fuel], "C1") || 
       !strcasecmp(Fuels[Fuel], "C2") ||
       !strcasecmp(Fuels[Fuel], "C3") ||
       !strcasecmp(Fuels[Fuel], "C4") ||
       !strcasecmp(Fuels[Fuel], "C5") ||
       !strcasecmp(Fuels[Fuel], "C7") ||
       !strcasecmp(Fuels[Fuel], "D1") ||
       !strcasecmp(Fuels[Fuel], "S1") ||  
       !strcasecmp(Fuels[Fuel], "S2") ||
       !strcasecmp(Fuels[Fuel], "S3")  )
      RSI = a[Fuel] * pow(1. - exp(-b[Fuel] * ISI), c[Fuel]);          /* 26 */

   if (!strcasecmp(Fuels[Fuel], "M1"))
      RSI =       PC/100. * ROScalc(FBPFuel("C2"), ISI, NoBUI, FMC, SFC, PC, PDF, C)
         + (100.-PC)/100. * ROScalc(FBPFuel("D1"), ISI, NoBUI, FMC, SFC, PC, PDF, C);  /* 27 */

   if (!strcasecmp(Fuels[Fuel], "M2"))
      RSI =          PC/100. * ROScalc(FBPFuel("C2"), ISI, NoBUI, FMC, SFC, PC, PDF, C)
      +0.2*(100.-PC)/100. * ROScalc(FBPFuel("D1"), ISI, NoBUI, FMC, SFC, PC, PDF, C);  /* 28 */

   if (!strcasecmp(Fuels[Fuel], "M3"))
   {
      a[Fuel] = 170. * exp(-35.0/PDF);                                 /* 29 */
      b[Fuel] = 0.082 *exp(-36.0/PDF);                                 /* 30 */
      c[Fuel] = 1.698 - 0.00303 * PDF;                                 /* 31 */
      RSI = a[Fuel] * pow(1. - exp(-b[Fuel] * ISI), c[Fuel]);          /* 26 */
   }

   if (!strcasecmp(Fuels[Fuel], "M4"))
   {
      a[Fuel] = 140. * exp(-33.5/PDF);      /* Changed from 35.5 to 33.5; 32 */
      b[Fuel] = 0.0404;                                                /* 33 */
      c[Fuel] = 3.02 * exp(-0.00714*PDF);                              /* 34 */
      RSI = a[Fuel] * pow(1. - exp(-b[Fuel] * ISI), c[Fuel]);          /* 26 */
   }

   if (!strcasecmp(Fuels[Fuel], "O1a") ||
       !strcasecmp(Fuels[Fuel], "O1b"))
   {
      if (C > 50.)
      {
         CF = 0.02 * C - 1.0;                                          /* 35 */
/* RSI has been substituted for ROS in eqn 36 */
         RSI = a[Fuel] * pow(1. - exp(-b[Fuel] * ISI), c[Fuel]) * CF;  /* 36 */
      }
      else
         RSI = 0.;
   }

   if (!strcasecmp(Fuels[Fuel], "C6"))
      ROS = C6calc(Fuel, ISI, BUI, FMC, SFC);
   else
   {
      BE = BEcalc(Fuel, BUI);
      ROS = BE*RSI;
   }

   if (ROS <= 0.) ROS = 0.000001;

   return(ROS);
}

double ISICalc(double FFMC, double WSV)

{
  double fW,m,fF,ISI;
     
     m = 147.2*(101.-FFMC)/(59.5+FFMC);                                /* 46 */
     fF = 91.9*exp(-0.1386*m)*(1.+pow(m,5.31)/4.93e7);                 /* 45 */

   if (WSV < 40.)
      fW = exp(0.05039*WSV);                                           /* 53 */
   else
      fW = 12. * (1. - exp(-0.0818 * (WSV-28.)));                     /* 53a */

     ISI = 0.208*fW*fF;                                                /* 52 */

  return(ISI);
}

/* Effect of Slope on Rate of Spread */

void Slopecalc(int Fuel, double FFMC, double BUI, double WS, double WAZ, double GS, 
               double SAZ, double FMC, double SFC, double PC, double PDF, double C,
               double *pRAZ, double *pWSV)
{
double SF, RSZ, RSF, ISZ, ISF, CF, m, fF, WSE, WSX, WSY; 

   if (GS > 60.) GS=60.;
   SF = exp(3.533 * pow(GS/100., 1.2));                                /* 39 */
   ISZ = ISICalc(FFMC, 0.);
   RSZ = ROScalc(Fuel, ISZ, NoBUI, FMC, SFC, PC, PDF, C);
   RSF = RSZ * SF;                                                     /* 40 */

   if (!strcasecmp(Fuels[Fuel], "C1") || 
       !strcasecmp(Fuels[Fuel], "C2") ||
       !strcasecmp(Fuels[Fuel], "C3") ||
       !strcasecmp(Fuels[Fuel], "C4") ||
       !strcasecmp(Fuels[Fuel], "C5") ||
       !strcasecmp(Fuels[Fuel], "C6") ||
       !strcasecmp(Fuels[Fuel], "C7") ||
       !strcasecmp(Fuels[Fuel], "D1") ||
       !strcasecmp(Fuels[Fuel], "S1") ||
       !strcasecmp(Fuels[Fuel], "S2") ||
       !strcasecmp(Fuels[Fuel], "S3"))
      ISF = log(1. - pow(RSF/a[Fuel], 1./c[Fuel]))/(-b[Fuel]);         /* 41 */

   if (!strcasecmp(Fuels[Fuel], "M1") ||
       !strcasecmp(Fuels[Fuel], "M2"))
      ISF = log(1. - pow(RSF/(PC/100.*a[FBPFuel("C2")])  , 1./c[FBPFuel("C2")] ) )
                    /(-b[FBPFuel("C2")] );                                 /* 42 */

   if (!strcasecmp(Fuels[Fuel], "M3"))
   {
      a[Fuel] = 170. * exp(-35.0/PDF);                                 /* 29 */
      b[Fuel] = 0.082 *exp(-36.0/PDF);                                 /* 30 */
      c[Fuel] = 1.698 - 0.00303 * PDF;                                 /* 31 */
      ISF = log(1. - pow(RSF/a[Fuel], 1./c[Fuel]))/(-b[Fuel]);         /* 41 */
   }

   if (!strcasecmp(Fuels[Fuel], "M4"))
   {
      a[Fuel] = 140. * exp(-33.5/PDF);                                 /* 32 */
      b[Fuel] = 0.0404;                                                /* 33 */
      c[Fuel] = 3.02 * exp(-0.00714*PDF);                              /* 34 */
      ISF = log(1. - pow(RSF/a[Fuel], 1./c[Fuel]))/(-b[Fuel]);         /* 41 */
   }

   if (!strcasecmp(Fuels[Fuel], "O1a") ||
       !strcasecmp(Fuels[Fuel], "O1b"))
       if (C > 50.)
       {
         CF = 0.02 * C - 1.0;                                          /* 35 */
         ISF = log(1. - pow(RSF/(CF*a[Fuel]), 1./c[Fuel]))/(-b[Fuel]); /* 43 */
       }
      else
         ISF = 0.;

   m = 147.2*(101.-FFMC)/(59.5+FFMC);                                  /* 46 */
   fF = 91.9*exp(-.1386*m)*(1.+pow(m,5.31)/4.93e7);                    /* 45 */

   WSE = log(ISF/(0.208 * fF))/0.05039;                                /* 44 */

   WSX = WS*sin(WAZ) + WSE*sin(SAZ);                                   /* 47 */
   WSY = WS*cos(WAZ) + WSE*cos(SAZ);                                   /* 48 */
   *pWSV = sqrt(WSX*WSX + WSY*WSY);                                    /* 49 */
   *pRAZ = acos(WSY/ (*pWSV));   /* in radians */                      /* 50 */
   if (WSX < 0.) *pRAZ = 2*PI - *pRAZ;                                 /* 51 */
}

/* BUI Effect */

double BEcalc(int Fuel, double BUI)
{
double BE;
   if (BUI > 0. && BUIo[Fuel] > 0.)
      BE = exp(50.*log(q[Fuel])*(1./BUI - 1./BUIo[Fuel]));             /* 54 */
   else
      BE = 1.;

   return (BE);
}

/* Crown Fraction Burned (CFB) calculation */

double CFBcalc(int Fuel, double FMC, double SFC, double ROS)
{
double CFB, CSI, RSO;

/* I have no idea why BMW has this in there 
   if (!strcasecmp(Fuels[Fuel], "C6")) SFC = 2.0;*/

   CFB = 0.;
   CSI = 0.001 * pow(CBH[Fuel], 1.5) * pow(460. + 25.9*FMC, 1.5);      /* 56 */
   RSO = CSI/(300.*SFC);                                               /* 57 */
   if (ROS > RSO) CFB = 1. - exp(-0.23*(ROS - RSO));                   /* 58 */
   return (CFB);
}

double C6calc(int Fuel, double ISI, double BUI, double FMC, double SFC)
{
double T, h, FME, RSI, RSS, RSC, FMEavg, CFB, ROS;

   FMEavg = 0.778;                                                /* page 37 */
   T = 1500. - 2.75 * FMC;                                             /* 59 */
   h = 460. + 25.9 * FMC;                                              /* 60 */
   FME = pow(1.5 - 0.00275 * FMC, 4.)/(460. + 25.9*FMC) * 1000.;       /* 61 */
   RSI = 30. * pow(1. - exp(-0.08 * ISI), 3.0);                        /* 62 */
   RSS = RSI * BEcalc(Fuel, BUI);                                      /* 63 */
   RSC = 60. * (1. - exp(-0.0497*ISI)) * FME/FMEavg;                   /* 64 */

   if (RSC > RSS)
   {
      CFB = CFBcalc(Fuel, FMC, SFC, RSS);
      ROS = RSS + CFB*(RSC-RSS);                                       /* 65 */
   }
   else
      ROS = RSS;

   return(ROS);
}

double TFCcalc(int Fuel, double CFB, double SFC, double PC)
{
double CFC, TFC;
   CFC = CFL[Fuel] * CFB;                                              /* 66 */
   if (!strcasecmp(Fuels[Fuel], "M1") ||
       !strcasecmp(Fuels[Fuel], "M2"))
      CFC = PC/100.*CFC;
   TFC = SFC + CFC;                                                    /* 67 */

   return(TFC);
}

double FIcalc(double FC, double ROS)
{
double FI;
   FI = 300. * FC * ROS;                                               /* 69 */

   return(FI);
}

double BROScalc (int Fuel, double FFMC, double BUI, double WSV, double FMC, double SFC, double PC, double PDF, double C)
{
double m, fF, BfW, BISI, BROS;

   m = 147.2*(101.-FFMC)/(59.5+FFMC);                                  /* 46 */
   fF = 91.9*exp(-0.1386*m)*(1.+pow(m,5.31)/4.93e7);                   /* 45 */
   BfW = exp(-0.05039*WSV);                                            /* 75 */
   BISI = 0.208*BfW*fF;                                                /* 76 */

/* Note the BUI effect is captured in ROScalc */
   BROS =  ROScalc(Fuel, BISI, BUI, FMC, SFC, PC, PDF, C);                      /* 77 */

   return (BROS);
}

double ROStcalc(int Fuel, double ROSeq, double t, double CFB)
{
double ROSt, alpha;

   if (!strcasecmp(Fuels[Fuel], "C1") ||
       !strcasecmp(Fuels[Fuel], "O1a") ||
       !strcasecmp(Fuels[Fuel], "O1b") ||
       !strcasecmp(Fuels[Fuel], "S1") ||
       !strcasecmp(Fuels[Fuel], "S2") ||
       !strcasecmp(Fuels[Fuel], "S3"))
      alpha = 0.115;                                              /* page 41 */
   else
      alpha = 0.115 - 18.8 * pow(CFB, 2.5) * exp(-8.* CFB);            /* 72 */

   ROSt = ROSeq * (1. - exp (-alpha * t));                             /* 70 */

   return(ROSt);
}

double LBcalc(int Fuel, double WSV)
{
double LB;

   if (!strcasecmp(Fuels[Fuel], "O1a") ||
       !strcasecmp(Fuels[Fuel], "O1b"))

      if (WSV >= 1.0)
         LB = 1.1 * pow(WSV, 0.464);        /* corrected from "+" to "*"; 80 */
      else
         LB = 1.0;                                                     /* 81 */

   else
      LB = 1.0 + 8.729 * pow(1.-exp(-0.030*WSV), 2.155);               /* 79 */

   return(LB);
}

double FROScalc(double ROS, double BROS, double LB)
{
double FROS;

   FROS = (ROS + BROS)/LB/2.;                                          /* 89 */

   return(FROS);
}

/* int Date2Julian(int month, int day, int year)
{
/* remember Jan=1, Feb=2, ... 
int days[12] = {0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334};
int julian;

   julian = days[month-1]+day;
   if (year%4 == 0  && year > 0 && month > 2) julian = julian+1;

   return(julian); 
} */

int FBPdefaults(FBP *fbp)
{
   strcpy(fbp->FuelType,"C2"); 
   fbp->Accel = 0;
   fbp->Dj = 0;
   fbp->Do = 0;
   fbp->ELV = -9;
   fbp->BUIEff = 1;
   fbp->t = 0;
   fbp->FFMC = 90;
   fbp->ISI = 0;
   fbp->BUI = 60;
   fbp->WS = 10;
   fbp->WD = 0;
   fbp->GS = 0;
   fbp->Aspect = 0;
   fbp->PC = 50;
   fbp->PDF = 35;
   fbp->C = 80.;
   fbp->GFL = 3.;
   fbp->CBH = 7.;
   fbp->LAT = 60.;
   fbp->LON = -120.;
   fbp->FMC = 0.;
   fbp->theta = 0.;

   return(1);
}

/*void FBPtest()
{
double SAZ;
int status, Date1, Date2, Do, month, day, year;

FBP fbp;

   infile  = fopen("fbp97.inp", "r");
   outfile = fopen("fbp97.out", "w");

   while (!feof(infile))
   {

      fscanf(infile, "%s%lf%d%d%d%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%d", 
             fbp.FuelType, &fbp.t, &fbp.Accel, &Date1, &Date2,
			 &fbp.LAT, &fbp.LON, &fbp.FFMC, &fbp.WS, &fbp.WD, &fbp.BUI, 
             &fbp.GS, &SAZ, &fbp.PC, &fbp.PDF, &fbp.GFL, &fbp.C, &fbp.ELV);

      printf("%s %lf %d %d %d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %d\n\n",  
             fbp.FuelType, fbp.t, fbp.Accel, Date1, Date2, 
			 fbp.LAT, fbp.LON, fbp.FFMC, fbp.WS, fbp.WD, fbp.BUI, 
             fbp.GS, SAZ, fbp.PC, fbp.PDF, fbp.GFL, fbp.C, fbp.ELV);

      if (fbp.Accel = 2) fbp.Accel = 0;
      fbp.LON = -fbp.LON;
      fbp.Aspect = SAZ-180.;
      fbp.BUIEff = 1;
      fbp.CBH = -1.;


      month = (int)(Date1/100);
      day = Date1 - 100*month;
      year = -1;
      fbp.Dj = Date2Julian(month, day, year);

      if (Date2 > 0)
      {
         month = (int)(Date2/100);
         day = Date2 - 100*month;
         year = -1;
         Do = Date2Julian(month, day, year);
      }
      else
         Do = -1;

      status = FBPCalc(&fbp);
      printf("%3s %5.2lf %3.1lf %5.1lf %5.2lf %6.0lf\n", fbp.FuelType, fbp.SFC, fbp.CFB, fbp.ROS, fbp.TFC, fbp.HFI);
      fprintf(outfile, "%3s %5.2lf %3.1lf %5.1lf %5.2lf %6.0lf\n", fbp.FuelType, fbp.SFC, fbp.CFB, fbp.ROS, fbp.TFC, fbp.HFI);

   }

}*/