#include <math.h>
#include <string.h>
#include <stdio.h>

/* SpotCalc calculates the maximum spot fire distance following

Albini, F.A.  1979.  Spot Fire Distance from Burning Trees - 
	A Predictive Model.  USDA For. Serv. Gen. Tech. Rep. INT-56, 
	Intermt. For. and Range Exp. Stn.  Ogden UT.  73 pp	

  and

Albini, F.A.  1983.  Potential Spotting Distance from Wind=Driven
	Surface Fires.  USDA For. Serv. Gen. Tech. Rep. INT-309, 
	Intermt. For. and Range Exp. Stn.  Ogden UT.  27 pp	
    
Chase, Carolyn H.  Spot fire distance equations for pocket calculators.  
    Res. Note INT-310, Intermt. For. and Range Exp. Stn.  Ogden UT.  19 pp	

*/

double SpotCalc (char * Fuel, double HFI, double WS, double H)

/*	Fuel is the FBP fuel type 
	HFI is the head fire intensity [kW/m]
	WS is the 10 metre wind speed [km/h]
	H is the tree height [m] */

{
double  X,  z, z_o, z_F, U, U_H, Uh, A, B, h, fU, E, t_o, t_T, a, b,
		h_c, h_star, F, S,	// these are from Chase '81
		g = 9.8,		// graviational accelleration
		z_10 = 10.,		// 10 m wind speed height
		D = 0.0254;	// 1 inch nominal firebrand diameter

int i, j, k;

//	printf("%s %lf %lf %lf\n", Fuel, HFI, WS, H);

	U = WS / 3.6;	/* convert to metres per second */ 

/* Spot fires from burning trees  - Albini 1979 */
	if (!strnicmp(Fuel, "C", 1) ||
		!strnicmp(Fuel, "M", 1) ||
		!strnicmp(Fuel, "D", 1) )

	{
		if (H <= 0.)
			H = 10.;	/* assume a default height of ten metres */

            
/* This represents the principle difference between Albini's model and this subroutine.
   By estimating the flame height from HFI, the program links Albini's work to the FBP system */

/* z_F = height of flame tip */
		z_F = sqrt(HFI/300.);												/* Byram (1959) */



/*** Appendix D: Lofting of Firebrands by Flame and Bouyant Plume ***/

		a = 5.963;		/* D33 */ 
		b = 4.563;		/* D34 */

/*  The following section was written to confirm the calculations match the results shown in Figure D-9

		for (i=0; i< 5; i++)
			for (j=0; j<8; j++)
			{
				z = (2 * i) * z_F;
				t_T = 2.* j;
				t_o = t_T - (1.2 + a/3.*(pow( (b+z/z_F)/a, 1.5)-1));		/* D43 * /
				printf("%d %d %lf\n", 2*i, 2*j, t_o);						/* Figure D-9 * /
			} 

		for (j=0; j<8; j++)
			for (k=4; k<12; k++)
			{
				z_o = 0.2*z_F;
				z = j*z_F;
				D = k*z_F/100000.;
				t_o = t_T - (1.2 + a/3.*(pow( (b+z/z_F)/a, 1.5)-1));		/* D43 * /

					printf("%d %d %d %lf %lf\n", i,j,k,t_o,t_T);			/* trying to reproduce D-10 * /
				}
*/

/* Now this is not clear...
   z, or z(0), = initial firebrand height above ground  (I believe represents the maximum height of the fire brand) 
   z_o = initial particle height, which is half the tree top height (based on table from page 53)? */
		z_o = 0.5 * H;
		z = z_F*( 3900 * (D/z_F) - z_o/z_F);								/* D45 */

/* This is better: Albini 1981 (page 4) for continuous flames */
		z = 12.2*z_F;

/*** Appendix E: A Crude Model for Surface Wind over Rought Terrain ***/

/* z_o = roughness height */
		z_o = 0.1313 * H;													/* E2 */
		U_H = U /0.493 / log(z_10/z_o);										/* E4 */

/*** Appendix F: Firebrand Trajectories ***/

/*	Fig 8a seems to work
		z=450.;
		H = 100.;
		g = 9.8 * 39.37/12;	// convert to ft/sec^2
		U_H = 25./3600; */	

/*	Fig 8b seems to work 
		z=180.;
		H = 40.;
		U_H = 40./3.6; */

/* z = z(0) */
		X = 21.9 * U_H * sqrt(H/g)* (0.362 + sqrt(z/H)*(0.5)* log(z/H) );	/* F22 */
	}

/* Spot fires from wind-driven surface fires - Albini 1983 */
	else
	{
						
		A=0.; B=0.; X = 0.;

		if (!strnicmp(Fuel, "S", 1))
		{
			A = 69.6; B = -0.818;	/* Medium conifer slash from Table 4 (Albini 1983) */
			h = 1.;					/* h is now the height of the vegetative cover */
		}

		if (!strnicmp(Fuel, "O", 1))
		{
			A = 129.; B = -1.238;	/* Tall grass from Table 4 */
			h = 0.5;				/* h is now the height of the vegetative cover */
		}

		if (!strnicmp(Fuel, "O1A", 3))
		{
			A = 162.; B = -1.465;	/* Tall grass from Table 4 */
			h = 1.;				/* h is now the height of the vegetative cover */

		}

		if (A > 0.)
		{
/*	Mature chaparral
    		A = 83.9;   
			B = -1.04;
			U = 20.;
			HFI = 50000.;*/


			fU = A*pow(U, B);												/* 8 */
//			printf("f(U) = %lf\n", fU);

			E = HFI* fU;													/* 7  */
//			printf("E = %lf\n", E);

			H = 0.173*sqrt(E);												/* 2 */
//			printf("H = %lf\n", H);

			h_c = pow(H, 0.337) - 1.22;										/* Chase 1981 */
//			printf("h_c = %lf\n", h_c);

			if (h_c < 6.) 
				h = 6.;
			else
				h = h_c;

			Uh = U*pow(h/10., 1./7);										/* 4 */
//			printf("Uh = %lf\n", 3.6*Uh);

			U_H = U*pow(H/10., 1./7);										/* 4 */
//			printf("U_H = %lf\n", U_H);

			h=h_c;

			F = 1.3 * 3.6 * Uh * sqrt(h) * (0.362 + sqrt(H/h) *(0.5)* log(H/h));	/* Chase '81 */
//			printf("\nF = %lf\n\n", F);

/* Downwind drift */
			if (h < 10.)
				X = 2.78 * U_H * sqrt(H);									/* 3 */
			else
				X = Uh * sqrt(H) * (2.13 + log (H/h));						/* 5 */

			X = X+F;
//			printf("X = %lf\n", X);

		}
	}
   return (X);
}