%Function ionocorr.m again from David Johnston, Essex Univ.,
%Modified by Keith Alter

function Iono = ionocorr(Latu, Lonu, Az, El, Ttr, alpha, beta)


        %/*for clarity copy array*/
        a0 = alpha(1);
        a1 = alpha(2);
        a2 = alpha(3);
        a3 = alpha(4);
        b0 = beta(1);
        b1 = beta(2);
        b2 = beta(3);
        b3 = beta(4);

        %/*convert from degrees to semicircles*/
        Latu = Latu / 180;
        Lonu = Lonu / 180;
	El = El / 180;
        Az = Az / 180;

	% Now "crunch" numbers per the Klobuchar model algorithm
	tsi = 0.0137/(El+0.11)-0.022;
	phi_l = Latu + tsi*cos(Az*pi);
	if phi_l >.416
		phi_l = .416;
	end
	if phi_l <-.416
		phi_l = -.416;
	end
	lambda_l = Lonu + (tsi*sin(Az*pi)/cos(phi_l*pi));
	phi_m = phi_l+0.064*cos((lambda_l-1.617)*pi);
	
	AMP = (a0+a1*phi_m+a2*phi_m^2+a3*phi_m^3);
	if AMP < 0
		AMP = 0;
	end

	PER = (b0+b1*phi_m+b2*phi_m^2+b3*phi_m^3);
	if PER < 72000
		PER = 72000;
	end

	t = 4.32e4*lambda_l+Ttr;
	t = mod(t,86400);
	F = 1.0+16*(0.53-El)^3;
	x = 2*pi*(t-50400)/PER;
	if abs(x) > 1.57
		Iono = F*5e-9;
	else
		Iono = F*(5e-9+AMP*(1-x^2/2+x^4/24));
	end;
	
