Roth1972.c Source File

00001 
00025 #include "Roth1972.h"
00026 
00027 #define ROTH_1972_WTG_CLASSES           (6)
00028 static const int   TimeLagClass         [ROTH_1972_WTG_CLASSES] = {0, 0, 1, 1, 2, 2};
00029 static const float WtgSzClassEnglish    [ROTH_1972_WTG_CLASSES] = {1200.0,  192.0,   96.0,   48.0,  16.0,  0.0};
00030 static const float WtgSzClassMetric     [ROTH_1972_WTG_CLASSES] = {3633.61, 581.37, 290.68, 145.34, 48.44, 0.0};
00031 
00032 int Roth1972FireSpreadSetFuelBed(RothFuelModel * rfm)   {
00033     double lload,   dload;
00034     double lhc,     dhc;
00035     double lseff,   dseff;
00036     double letas,   detas;
00037 
00038     double flive, beta_opt, ratio, aa, sigma_15, gamma, gamma_max, c, e;
00039 
00040     /* check args */
00041     if ( rfm == NULL )  {
00042         ERR_ERROR("RothFuelModel not initialized, FuelBed step in FireSpread Pipeline failed. \n", ERR_EINVAL);
00043         }
00044     if ( rfm->rp == NULL )  {
00045         ERR_ERROR("RothPipeline not initialized, FuelBed step in FireSpread Pipeline failed. \n", ERR_EINVAL);
00046         }
00047     
00048     /* return if RothFuelModel represents unburnable fuel */    
00049     if ( rfm->brntype == EnumRothUnBurnable )   {
00050         return ERR_SUCCESS;
00051         }
00052     
00053     /* reset pipeline to defualt state */
00054     if ( RothPipelineSetDefaultValues(rfm->rp) )    {
00055         ERR_ERROR("FuelBed step in FireSpread Pipeline failed. \n", ERR_EBADFUNC);
00056         }
00057     
00058     /* set pipeline state */
00059     rfm->rp->pipe = EnumSetFuelBedPipe;
00060         
00061     /* all calculations done in ENGLISH units */
00062     if ( rfm->units == EnumMetricUnits )    {
00063         if ( RothFuelModelMetricToEnglish(rfm) )    {
00064             ERR_ERROR("FuelBed step in FireSpread Pipeline failed. \n", ERR_EBADFUNC);
00065             }
00066         }
00067     
00068     /* if fuel bed has no surface area go no further */
00069     if ( UNITS_FP_IS_ZERO(rfm->larea + rfm->darea) )    {
00070         return ERR_SUCCESS;
00071         }
00072 
00073     /* initialize scratch vars */
00074     flive = beta_opt = ratio = aa = sigma_15 = gamma = gamma_max = c = e = 0.0;
00075             
00076     /* apply weights to generate fuel bed properties */
00077     lload = (rfm->awtg[EnumLHSizeClass] * rfm->load[EnumLHSizeClass] * (1.0 - rfm->stot[EnumLHSizeClass])) +
00078             (rfm->awtg[EnumLWSizeClass] * rfm->load[EnumLWSizeClass] * (1.0 - rfm->stot[EnumLWSizeClass]));
00079     dload = (rfm->awtg[EnumD1HSizeClass] * rfm->load[EnumD1HSizeClass] * (1.0 - rfm->stot[EnumD1HSizeClass])) +
00080             (rfm->awtg[EnumD10HSizeClass] * rfm->load[EnumD10HSizeClass] * (1.0 - rfm->stot[EnumD10HSizeClass])) +
00081             (rfm->awtg[EnumD100HSizeClass] * rfm->load[EnumD100HSizeClass] * (1.0 - rfm->stot[EnumD100HSizeClass]));
00082             
00083     lhc =   (rfm->awtg[EnumLHSizeClass] * rfm->hc[EnumLHSizeClass]) +
00084             (rfm->awtg[EnumLWSizeClass] * rfm->hc[EnumLWSizeClass]);
00085     dhc =   (rfm->awtg[EnumD1HSizeClass] * rfm->hc[EnumD1HSizeClass]) +
00086             (rfm->awtg[EnumD10HSizeClass] * rfm->hc[EnumD10HSizeClass]) +
00087             (rfm->awtg[EnumD100HSizeClass] * rfm->hc[EnumD100HSizeClass]);
00088             
00089     lseff = (rfm->awtg[EnumLHSizeClass] * rfm->seff[EnumLHSizeClass]) +
00090             (rfm->awtg[EnumLWSizeClass] * rfm->seff[EnumLWSizeClass]);
00091     dseff = (rfm->awtg[EnumD1HSizeClass] * rfm->seff[EnumD1HSizeClass]) +
00092             (rfm->awtg[EnumD10HSizeClass] * rfm->seff[EnumD10HSizeClass]) +
00093             (rfm->awtg[EnumD100HSizeClass] * rfm->seff[EnumD100HSizeClass]);
00094 
00095     /* etas */
00096     letas = 1.0;
00097     if ( UNITS_FP_GT_ZERO(lseff) )  {
00098         letas = 0.174 / pow(lseff, 0.19);
00099         if ( letas > 1.0 )  {
00100             letas = 1.0;
00101             }
00102         }
00103     detas = 1.0;
00104     if ( UNITS_FP_GT_ZERO(dseff) )  {
00105         detas = 0.174 / pow(dseff, 0.19);
00106         if ( detas > 1.0 )  {
00107             detas = 1.0;
00108             }
00109         }
00110         
00111     /* rx factors */
00112     rfm->rp->lrx = lload * lhc * letas; 
00113     rfm->rp->drx = dload * dhc * detas;
00114     
00115     /* residence time */
00116     rfm->rp->taur = 384.0 / rfm->fsav;
00117     
00118     /* propagating flux */
00119     rfm->rp->ppflux = exp((0.792 + 0.681 * sqrt(rfm->fsav)) * (rfm->pkrat + 0.1)) / (192.0 + 0.2595 * rfm->fsav);
00120     
00121     /* gamma */
00122     beta_opt = 3.348 / (pow(rfm->fsav, 0.8189));
00123     ratio = rfm->pkrat / beta_opt;
00124     aa = 133.0 / (pow(rfm->fsav, 0.7913));
00125     sigma_15 = pow(rfm->fsav, 1.5);
00126     gamma_max = sigma_15 / (495.0 + 0.0594 * sigma_15);
00127     gamma = gamma_max * pow(ratio, aa) * exp(aa * (1.0 - ratio));   
00128 
00129     /* factor gamma into rx factors */
00130     rfm->rp->lrx *= gamma;
00131     rfm->rp->drx *= gamma;
00132     
00133     /* slope and wind intermediaries */
00134     rfm->rp->slp_k = 5.275 * pow(rfm->pkrat, -0.3);
00135     rfm->rp->wnd_b = 0.02526 * pow(rfm->fsav, 0.54);
00136     c = 7.47 * exp(-0.133 * pow(rfm->fsav, 0.55));
00137     e = 0.715 * exp(-0.000359 * rfm->fsav);
00138     rfm->rp->wnd_k = c * pow(ratio, -e);
00139     rfm->rp->wnd_e = pow(ratio, e) / c; 
00140     
00141     /* if no live fuels in bed go no further */
00142     if ( UNITS_FP_IS_ZERO(lload) )  {
00143         return ERR_SUCCESS;
00144         }
00145         
00146     /* live and dead fuel ratio and factor */
00147     if ( UNITS_FP_GT_ZERO(rfm->sav[EnumLHSizeClass]) )  {
00148         flive += rfm->load[EnumLHSizeClass] * exp(-500.0 / rfm->sav[EnumLHSizeClass]);
00149         }
00150     if ( UNITS_FP_GT_ZERO(rfm->sav[EnumLWSizeClass]) )  {
00151         flive += rfm->load[EnumLWSizeClass] * exp(-500.0 / rfm->sav[EnumLWSizeClass]);
00152         }
00153     rfm->rp->fdead = (rfm->load[EnumD1HSizeClass]  * rfm->effhn[EnumD1HSizeClass]) +
00154                      (rfm->load[EnumD10HSizeClass] * rfm->effhn[EnumD10HSizeClass]) +
00155                      (rfm->load[EnumD100HSizeClass] * rfm->effhn[EnumD100HSizeClass]);
00156     if ( UNITS_FP_GT_ZERO(flive) )  {
00157         rfm->rp->lmex = 2.9 * rfm->rp->fdead / flive;
00158         }                    
00159                 
00160     return ERR_SUCCESS;
00161     }
00162 
00163 int Roth1972FireSpreadNoWindNoSlope(RothFuelModel * rfm, double d1hfm, double d10hfm, double d100hfm,
00164                                                         double lhfm, double lwfm)   {
00165     double letam,   detam;
00166     double lm,      dm;
00167     double lmex,    dmex;
00168                                             
00169     double wfmd, rbqig, fdmois, qig, ratio;
00170     
00171     int     sz_cls  [EnumNumSizeClasses] = {0};
00172     double  tlag_cls[EnumNumSizeClasses] = {0.0};
00173     int i, j;
00174     
00175     /* check args */
00176     if ( rfm == NULL )  {
00177         ERR_ERROR("RothFuelModel not initialized, NoWindNoSlope step in FireSpread Pipeline failed. \n", ERR_EINVAL);
00178         }
00179     if ( rfm->rp == NULL )  {
00180         ERR_ERROR("RothPipeline not initialized, NoWindNoSlope step in FireSpread Pipeline failed. \n", ERR_EINVAL);
00181         }
00182     
00183     /* return if RothFuelModel represents unburnable fuel */    
00184     if ( rfm->brntype == EnumRothUnBurnable )   {
00185         return ERR_SUCCESS;
00186         }
00187     
00188     /* check and set pipeline state */
00189     if ( rfm->rp->pipe < EnumSetFuelBedPipe )   {
00190         ERR_ERROR("SetFuelBed step not complete, unable to execute NoWindNoSlope. \n", ERR_ESANITY);
00191         }
00192     rfm->rp->pipe = EnumNoWindNoSlopePipe;
00193     
00194     /* check for change in moisture */
00195     if ( (UNITS_FP_ARE_EQUAL(rfm->rp->d1hfm, d1hfm)) && (UNITS_FP_ARE_EQUAL(rfm->rp->d10hfm, d10hfm)) &&
00196          (UNITS_FP_ARE_EQUAL(rfm->rp->d100hfm, d100hfm)) && (UNITS_FP_ARE_EQUAL(rfm->rp->lhfm, lhfm)) &&
00197          (UNITS_FP_ARE_EQUAL(rfm->rp->lwfm, lwfm)) )    {
00198         return ERR_SUCCESS;
00199         }
00200     
00201     /* assign new moistures */
00202     rfm->rp->d1hfm = d1hfm;
00203     rfm->rp->d10hfm = d10hfm;
00204     rfm->rp->d100hfm = d100hfm;
00205     rfm->rp->lhfm = lhfm;
00206     rfm->rp->lwfm = lwfm;
00207     
00208     /* zero any previously set pipeline vars */
00209     rfm->rp->ros_0 = rfm->rp->hpua = rfm->rp->rxint = 0.0;
00210     rfm->rp->ros_max = rfm->rp->ros_az_max = 0.0;
00211     rfm->rp->ros_any = rfm->rp->ros_az_any = 0.0;
00212 
00213     /* initialize scratch vars */
00214     wfmd = rbqig = fdmois = qig = ratio = 0.0;
00215     letam = detam = lm = dm = lmex = dmex = 0.0;
00216 
00217     /* accumulate moisture based upon size classes */
00218     for(i = 0; i < EnumNumSizeClasses; i++)     {
00219         for(j = 0; j < ROTH_1972_WTG_CLASSES; j++)  {
00220             if ( rfm->units == EnumEnglishUnits && rfm->sav[i] > WtgSzClassEnglish[j] ) {
00221                 sz_cls[i] = j;
00222                 break;
00223                 }
00224             if ( rfm->units == EnumMetricUnits && rfm->sav[i] > WtgSzClassMetric[j] )   {
00225                 sz_cls[i] = j;          
00226                 break;
00227                 }
00228             }
00229         }
00230     for(i = 0; i < EnumNumSizeClasses; i++)     {
00231         switch(TimeLagClass[sz_cls[i]]) {
00232             case 0:
00233                 tlag_cls[i] = d1hfm;
00234                 break;
00235             case 1:
00236                 tlag_cls[i] = d10hfm;
00237                 break;
00238             case 2:
00239                 tlag_cls[i] = d100hfm;
00240                 break;
00241             default:
00242                 ERR_ERROR_CONTINUE("Unrecognized TimeLagClass, NoWindNoSlope in FireSpread Pipeline. \n", ERR_WARNING);
00243                 tlag_cls[i] = 0.0;
00244                 break;
00245             }       
00246         }
00247     wfmd =  (tlag_cls[EnumD1HSizeClass] * rfm->effhn[EnumD1HSizeClass] * rfm->load[EnumD1HSizeClass]) +
00248             (tlag_cls[EnumD10HSizeClass] * rfm->effhn[EnumD10HSizeClass] * rfm->load[EnumD10HSizeClass]) +
00249             (tlag_cls[EnumD100HSizeClass] * rfm->effhn[EnumD100HSizeClass] * rfm->load[EnumD100HSizeClass]);
00250     rfm->fm[EnumD1HSizeClass] = tlag_cls[EnumD1HSizeClass];
00251     rfm->fm[EnumD10HSizeClass] = tlag_cls[EnumD10HSizeClass];
00252     rfm->fm[EnumD100HSizeClass] = tlag_cls[EnumD100HSizeClass]; 
00253     rfm->fm[EnumLHSizeClass] = lhfm;
00254     rfm->fm[EnumLWSizeClass] = lwfm;
00255     
00256     /* live fuel extinction moisture */
00257     if ( UNITS_FP_GT_ZERO(rfm->load[EnumLHSizeClass]) || UNITS_FP_GT_ZERO(rfm->load[EnumLWSizeClass]) ) {
00258         if ( UNITS_FP_GT_ZERO(rfm->rp->fdead) )     {
00259             fdmois = wfmd / rfm->rp->fdead;
00260             }
00261         else    {
00262             fdmois = 0.0;
00263             }
00264         lmex = ((rfm->rp->lmex * (1.0 - fdmois / rfm->mex)) - 0.226);
00265         if ( lmex < rfm->mex )  {
00266             lmex = rfm->mex;
00267             }
00268         }
00269         
00270     /* dead fuel extinction moisture */     
00271     dmex = rfm->mex;
00272 
00273     /* accumulate category weighted moisture */
00274     for(i = 0; i < EnumNumSizeClasses; i++) {
00275         qig = 250.0 + 1116.0 * rfm->fm[i];
00276         if ( i == EnumLHSizeClass || i == EnumLWSizeClass ) {
00277             rbqig += qig * rfm->awtg[i] * rfm->larea * rfm->effhn[i];
00278             }
00279         else    {
00280             rbqig += qig * rfm->awtg[i] * rfm->darea * rfm->effhn[i];
00281             }
00282         }
00283     rbqig *= rfm->fdens;
00284     lm = (rfm->awtg[EnumLHSizeClass] * rfm->fm[EnumLHSizeClass]) +
00285          (rfm->awtg[EnumLWSizeClass] * rfm->fm[EnumLWSizeClass]);
00286     dm = (rfm->awtg[EnumD1HSizeClass] * rfm->fm[EnumD1HSizeClass]) +
00287          (rfm->awtg[EnumD10HSizeClass] * rfm->fm[EnumD10HSizeClass]) +
00288          (rfm->awtg[EnumD100HSizeClass] * rfm->fm[EnumD100HSizeClass]);
00289 
00290     /* reaction intensity contributed by live fuels */
00291     ratio = 0.0;
00292     if ( UNITS_FP_GT_ZERO(lmex) )   {
00293         ratio = lm / lmex;
00294         letam = 1.0 - 2.59 * ratio + 5.11 * ratio * ratio - 3.52 * ratio * ratio * ratio;
00295         }
00296     if ( lm >= lmex )   {
00297         letam = 0.0;
00298         }
00299     rfm->rp->rxint += rfm->rp->lrx * letam;
00300     
00301     /* reaction intensity contributed by dead fuels */
00302     ratio = 0.0;
00303     if ( UNITS_FP_GT_ZERO(dmex) )   {
00304         ratio = dm / dmex;
00305         detam = 1.0 - 2.59 * ratio + 5.11 * ratio * ratio - 3.52 * ratio * ratio * ratio;
00306         }
00307     if ( dm >= dmex )   {
00308         detam = 0.0;
00309         }
00310     rfm->rp->rxint += rfm->rp->drx * detam;
00311     
00312     /* heat per unit area */
00313     rfm->rp->hpua = rfm->rp->rxint * rfm->rp->taur;
00314     
00315     /* no wind no slope ros */
00316     if ( UNITS_FP_GT_ZERO(rbqig) )  {
00317         rfm->rp->ros_0 = rfm->rp->rxint * rfm->rp->ppflux / rbqig;
00318         }
00319     else    {
00320         rfm->rp->ros_0 = 0.0;
00321         }
00322         
00323     /* zero other ros vars */
00324     rfm->rp->ros_max = rfm->rp->ros_0;
00325     rfm->rp->ros_any = rfm->rp->ros_0;
00326     rfm->rp->ros_az_max = rfm->rp->ros_az_any = 0.0;
00327     
00328     return ERR_SUCCESS;                                 
00329     }
00330     
00331 int Roth1972FireSpreadWindSlopeMax(RothFuelModel * rfm, double wnd_fpm, double wnd_az, 
00332                                                         double slp_pcnt, double asp)    {
00333     double upslp, az_max, phi_ew;
00334     double split_deg, split_rad;
00335     double slp_rate, wnd_rate, rv, spread_max;
00336     double x, y, al, a;
00337     double max_wnd, eff_wnd, lw_ratio, eccen;
00338     int do_eff_wnd, ck_wnd_lim, wnd_lim;
00339     
00340     /* check args */
00341     if ( rfm == NULL )  {
00342         ERR_ERROR("RothFuelModel not initialized, WindSlopeMax step in FireSpread Pipeline failed. \n", ERR_EINVAL);
00343         }
00344     if ( rfm->rp == NULL )  {
00345         ERR_ERROR("RothPipeline not initialized, WindSlopeMax step in FireSpread Pipeline failed. \n", ERR_EINVAL);
00346         }
00347     
00348     /* return if RothFuelModel represents unburnable fuel */    
00349     if ( rfm->brntype == EnumRothUnBurnable )   {
00350         return ERR_SUCCESS;
00351         }
00352     
00353     /* check and set pipeline state */
00354     if ( rfm->rp->pipe < EnumNoWindNoSlopePipe )    {
00355         ERR_ERROR("NoWindNoSlope step not complete, unable to execute WindSlopeMax. \n", ERR_ESANITY);
00356         }
00357     rfm->rp->pipe = EnumWindSlopeMaxPipe;
00358 
00359     /* convert slope to rise/run and check for change in slope */
00360     if ( UNITS_FP_LT_ZERO(slp_pcnt) )   {
00361         slp_pcnt = 0.0;
00362         }
00363     slp_pcnt /= 100.0;
00364     if ( !UNITS_FP_ARE_EQUAL(rfm->rp->slp, slp_pcnt) )  {
00365         rfm->rp->phi_s = rfm->rp->slp_k * slp_pcnt * slp_pcnt;
00366         rfm->rp->slp = slp_pcnt;
00367         }
00368         
00369     /* convert wind direction from 'out of' to 'to' */
00370     wnd_az = ((int)(wnd_az + 180.0)) % 360;
00371     
00372     /* check for change in windspeed */
00373     if ( !UNITS_FP_ARE_EQUAL(rfm->rp->wnd_fpm, wnd_fpm) )   {
00374         if ( UNITS_FP_GT_ZERO(wnd_fpm) )    {
00375             rfm->rp->phi_w = rfm->rp->wnd_k * pow(wnd_fpm, rfm->rp->wnd_b);
00376             }
00377         else    {
00378             rfm->rp->phi_w = 0.0;
00379             }
00380         rfm->rp->wnd_fpm = wnd_fpm;
00381         }
00382         
00383     /* combine wind and slope */
00384     phi_ew = rfm->rp->phi_s + rfm->rp->phi_w;
00385     wnd_lim = 0;
00386     lw_ratio = 1.0;
00387     eccen = 0.0;
00388     if (asp >= 180.0)   {
00389         upslp = asp - 180.0;
00390         }
00391     else    {
00392         upslp = asp + 180.0;
00393         }
00394         
00395     /* Situation 1: no fire spread or reaction intensity */
00396     if ( !UNITS_FP_GT_ZERO(rfm->rp->ros_0) )        {
00397         spread_max = az_max = 0.0;
00398         eff_wnd = 0.0;
00399         do_eff_wnd = ck_wnd_lim = 0;
00400         }
00401     /* Situation 2: no wind and no wind */
00402     else if ( !UNITS_FP_GT_ZERO(phi_ew) )   {
00403         phi_ew = eff_wnd = az_max = 0.0;
00404         spread_max = rfm->rp->ros_0;
00405         do_eff_wnd = ck_wnd_lim = 0;
00406         }
00407     /* Situation 3: wind with no slope */
00408     else if ( !UNITS_FP_GT_ZERO(slp_pcnt) ) {        
00409         eff_wnd = wnd_fpm;
00410         do_eff_wnd = 0;
00411         spread_max = rfm->rp->ros_0 * (1.0 + phi_ew);
00412         az_max = wnd_az;
00413         ck_wnd_lim = 1;
00414         }
00415     /* Situation 4: slope with no wind */
00416     else if ( !UNITS_FP_GT_ZERO(wnd_fpm) )  {
00417         spread_max = rfm->rp->ros_0 * (1.0 + phi_ew);
00418         az_max = upslp;
00419         do_eff_wnd = ck_wnd_lim = 1;
00420         }
00421     /* Situation 5: wind blows upslope */
00422     else if ( UNITS_FP_ARE_EQUAL(upslp, wnd_az) )   {        
00423         spread_max = rfm->rp->ros_0 * (1.0 + phi_ew);
00424         az_max = upslp;        
00425         do_eff_wnd = ck_wnd_lim = 1;        
00426         }
00427     /* Situation 6: wind blows cross slope */
00428     else    {
00429         /* recalculate spread rate in the optimal direction */
00430         if ( upslp <= wnd_az )  {
00431             split_deg = wnd_az - upslp;
00432             }
00433         else    {
00434             split_deg = 360.0 - upslp + wnd_az;
00435             }
00436         split_rad = ROTH_1972_DEG_TO_RAD(split_deg);
00437         slp_rate = rfm->rp->ros_0 * rfm->rp->phi_s;
00438         wnd_rate = rfm->rp->ros_0 * rfm->rp->phi_w;
00439         x = slp_rate + wnd_rate * cos(split_rad);
00440         y = wnd_rate * sin(split_rad);
00441         rv = sqrt( (x * x) + (y * y));
00442         spread_max = rfm->rp->ros_0 + rv;
00443 
00444         /* recalculate phi_ew in the optimal direction */
00445         phi_ew = spread_max / rfm->rp->ros_0 - 1.0;
00446         if ( UNITS_FP_GT_ZERO(phi_ew) ) {
00447             do_eff_wnd = 1;
00448             }
00449         else    {
00450             do_eff_wnd = 0;
00451             }
00452         ck_wnd_lim = 1;
00453 
00454         /* recalculate direction of maximum spread in azimuth degrees */
00455         al = asin(fabs(y) / rv);
00456         if ( x >= 0.0 ) {
00457             if ( y >= 0.0 ) {
00458                 a = al;
00459                 }
00460             else    {
00461                 a = ROTH_1972_PI + ROTH_1972_PI - al;
00462                 }
00463             }
00464         else    {
00465             if ( y >= 0.0 ) {
00466                 a = ROTH_1972_PI - al;
00467                 }
00468             else    {
00469                 a = ROTH_1972_PI + al;
00470                 }
00471             }
00472         split_deg = ROTH_1972_RAD_TO_DEG(a);
00473         az_max = upslp + split_deg;
00474         if ( az_max > 360.0 )   {
00475             az_max -= 360.0;
00476             }
00477         }
00478         
00479     /* recalculate effective windspeed based upon phi_ew */
00480     if ( do_eff_wnd )   {
00481         eff_wnd = pow((phi_ew * rfm->rp->wnd_e), (1.0 / rfm->rp->wnd_b));
00482         }
00483         
00484     /* if effective windspeed excedes max windspeed, scale back */
00485     if ( ck_wnd_lim )   {
00486         max_wnd = 0.9 * rfm->rp->rxint;
00487         if ( eff_wnd > max_wnd )    {
00488             if ( !UNITS_FP_GT_ZERO(max_wnd) )   {
00489                 phi_ew = 0.0;
00490                 }
00491             else    {
00492                 phi_ew = rfm->rp->wnd_k * pow(max_wnd, rfm->rp->wnd_b);
00493                 }
00494             spread_max = rfm->rp->ros_0 * (1.0 + phi_ew);
00495             eff_wnd = max_wnd;
00496             wnd_lim = 1;
00497             }
00498         }
00499     
00500     /* determine fire ellipse parameters */
00501     if ( UNITS_FP_GT_ZERO(eff_wnd) )    {
00502         lw_ratio = 1.0 + 0.002840909 * eff_wnd;
00503         eccen = sqrt(lw_ratio * lw_ratio - 1.0) / lw_ratio;
00504         }
00505         
00506     /* store results */
00507     rfm->rp->asp = asp;
00508     rfm->rp->wnd_vec = wnd_az;
00509     rfm->rp->phi_ew = phi_ew;
00510     rfm->rp->wnd_eff = eff_wnd;
00511     rfm->rp->wnd_lim = wnd_lim;
00512     rfm->rp->ros_max = rfm->rp->ros_any = spread_max;
00513     rfm->rp->ros_az_max = rfm->rp->ros_az_any = az_max;
00514     rfm->rp->lwratio = lw_ratio;
00515     rfm->rp->eccen = eccen;
00516         
00517     return ERR_SUCCESS;
00518     }
00519 
00520 int Roth1972FireSpreadGetAtAzimuth(RothFuelModel * rfm, double az)  {
00521     double dir_deg, dir_rad;
00522     
00523     /* check args */
00524     if ( rfm == NULL )  {
00525         ERR_ERROR("RothFuelModel not initialized, GetAtAzimuth step in FireSpread Pipeline failed. \n", ERR_EINVAL);
00526         }
00527     if ( rfm->rp == NULL )  {
00528         ERR_ERROR("RothPipeline not initialized, GetAtAzimuth step in FireSpread Pipeline failed. \n", ERR_EINVAL);
00529         }
00530     
00531     /* return if RothFuelModel represents unburnable fuel */    
00532     if ( rfm->brntype == EnumRothUnBurnable )   {
00533         return ERR_SUCCESS;
00534         }
00535     
00536     /* check and set pipeline state */
00537     if ( rfm->rp->pipe < EnumWindSlopeMaxPipe ) {
00538         ERR_ERROR("WindSlopeMax step not complete, unable to execute GetAtAzimuth. \n", ERR_ESANITY);
00539         }
00540     rfm->rp->pipe = EnumGetAtAzimuthPipe;
00541 
00542     /* no fire spread */
00543     if ( !UNITS_FP_GT_ZERO(rfm->rp->ros_max) )  {
00544         return ERR_SUCCESS;
00545         }
00546         
00547     /* combined wind slope factor is 0 or azimuth in direction of max spread */
00548     if ( !UNITS_FP_GT_ZERO(rfm->rp->phi_ew) || UNITS_FP_ARE_EQUAL(rfm->rp->ros_az_max, az) )    {
00549         rfm->rp->ros_any = rfm->rp->ros_max;
00550         }
00551     /* azimuth not in direction of max spread */
00552     else    {
00553         /* calculate angle btwn max ros and requested az */
00554         if ( (dir_deg = fabs(rfm->rp->ros_az_max - az)) > 180.0 )   {
00555             dir_deg = 360.0 - dir_deg;
00556             }
00557         dir_rad = ROTH_1972_DEG_TO_RAD(dir_deg);
00558         /* calculate ros in this direction */
00559         rfm->rp->ros_any = rfm->rp->ros_max * (1.0 - rfm->rp->eccen) / (1.0 - rfm->rp->eccen * cos(dir_rad));
00560         }
00561 
00562     /* store azimuth */     
00563     rfm->rp->ros_az_any = az;
00564     
00565     return ERR_SUCCESS;
00566     }
00567  
00568 /* end of Roth1972.c */