function [ cgmlat,ab_mlat ] = AURORAL_BOUNDARY( xkp,xmlt )
%AURORAL_BOUNDARY equatorward auroral boundary
%		subroutine auroral_boundary(xkp,xmlt,cgmlat,ab_mlat)
%-----------------------------------------------------------------------
% Computes equatorward auroral boundary values for givern kp value.
% kp given in units of 0.1 (xkp) for the range from 0.0 to 9.0. Model 
% values are only used for kp=0,1,2,3,4,5,6,7,8,9 and a linear inter-
% polation is applied for intermediate kp values.
% 
% The auroral oval boundary is given as an array for corrected magnetic 
% latitude CGM (ab_mlat). The 48 values correspond to the MLT values 
% of 0.0,0.5,1.0,1.5,2.0 .. 23.5. If the input xmlt is greater than
% -1 then the program determines the CGM latitude, cgmlat, that 
% corresponds to the given MLT value (xmlt).
%
% Y. Zhang and L.J. Paxton, An empirical Kp-dependent global auroral 
% model based on TIMED/GUVI FUV data, Journal of Atmospheric and 
% Solar-Terrestrial Physics 70, 1231–1242, 2008.
% 
%----------------------------------------------------------------------- 

%        dimension zp_mlat(48,10),ab_mlat(48),ab_mlt(48)
  persistent zp_mlat ab_mlt NI NJ;
  if isempty(zp_mlat)
    dat = [ ...
       66.1,65.9,65.9,66.0,66.2,66.5,66.7,67.0,67.4,67.8,68.2,68.7, ...
       69.4,70.0,70.5,70.8,71.0,71.3,71.9,72.7,73.8,75.4,77.0,77.8, ...
       77.3,76.7,76.4,76.2,75.9,75.4,74.7,74.1,73.5,73.0,72.5,71.9, ...
       71.6,71.1,70.6,70.0,69.2,68.5,67.9,67.5,67.2,66.9,66.7,66.4, ...
       63.0,62.9,63.0,63.0,63.1,63.2,63.3,63.6,64.0,64.5,65.0,65.6, ...
       66.2,66.7,67.1,67.6,68.1,68.6,69.4,70.3,71.5,72.8,74.1,74.9, ...
       74.8,74.7,74.7,74.5,74.0,73.1,72.2,71.2,70.4,69.7,68.9,68.2, ...
       67.5,66.9,66.3,65.8,65.4,65.0,64.5,64.1,63.7,63.3,63.1,62.9, ...
       61.1,61.3,61.5,61.7,61.9,62.1,62.2,62.4,62.7,63.0,63.4,64.0, ...
       64.4,65.0,65.5,65.9,66.4,66.9,67.5,68.3,69.3,70.5,71.6,72.3, ...
       72.7,72.8,72.8,72.4,71.6,70.6,69.7,68.9,68.2,67.5,66.7,65.9, ...
       65.0,64.3,63.7,63.4,63.1,62.8,62.4,62.0,61.5,61.2,61.0,61.0, ...
       59.6,60.0,60.4,60.7,60.7,60.7,60.5,60.4,60.6,61.1,61.8,62.5, ...
       63.0,63.4,63.7,64.1,64.6,65.4,66.2,67.0,67.7,68.5,69.3,70.1, ...
       70.6,70.9,70.9,70.4,69.3,68.0,66.9,66.0,65.2,64.5,63.7,63.0, ...
       62.3,61.7,61.2,60.9,60.6,60.4,60.2,60.1,59.8,59.6,59.4,59.4, ...
       58.5,58.8,59.2,59.4,59.4,59.2,58.9,58.9,59.2,59.7,60.5,61.2, ...
       61.7,62.0,62.3,62.6,63.3,64.1,65.1,65.9,66.4,66.9,67.5,68.4, ...
       69.0,69.2,68.9,68.1,66.7,65.3,64.2,63.4,62.7,62.0,61.1,60.4, ...
       59.7,59.2,58.9,58.6,58.4,58.4,58.4,58.5,58.5,58.4,58.3,58.3, ...
       57.6,57.8,57.9,57.9,57.7,57.5,57.4,57.5,57.9,58.6,59.3,59.8, ...
       60.3,60.5,60.8,61.2,62.0,63.0,64.1,64.9,65.4,65.7,66.2,66.8, ...
       67.2,66.8,66.0,64.8,63.4,62.1,61.3,60.7,60.1,59.3,58.2,57.1, ...
       56.6,56.3,56.1,56.0,56.0,56.2,56.6,56.9,57.1,57.3,57.4,57.5, ...
       54.3,54.9,55.4,55.6,55.6,55.3,55.1,55.0,55.3,55.8,56.5,57.2, ...
       57.8,58.3,58.7,59.4,60.3,61.2,62.0,62.8,63.4,63.9,64.1,64.1, ...
       63.8,63.3,62.5,61.5,60.2,58.9,57.7,56.7,56.0,55.5,55.1,54.8, ...
       54.6,54.3,53.9,53.4,53.2,53.1,53.3,53.4,53.5,53.5,53.6,53.8, ...
       52.9,53.6,54.2,54.6,54.6,54.3,54.0,53.7,53.8,54.1,54.8,55.7, ...
       56.4,56.9,57.4,58.0,58.8,59.6,60.2,60.8,61.7,62.4,62.6,62.2, ...
       61.5,60.7,59.8,58.9,57.9,56.8,55.5,54.5,53.8,53.4,53.3,53.3, ...
       53.5,53.2,52.7,52.1,51.7,51.7,51.7,51.9,51.9,52.0,52.1,52.4, ...
       51.8,52.5,53.2,53.6,53.6,53.3,52.9,52.5,52.4,52.6,53.3,54.2, ...
       55.0,55.6,56.1,56.7,57.3,57.9,58.1,58.6,59.7,60.7,60.9,60.3, ...
       59.2,58.1,57.1,56.3,55.5,54.6,53.5,52.5,51.8,51.4,51.4,51.7, ...
       52.0,51.9,51.4,50.8,50.4,50.3,50.4,50.5,50.6,50.7,50.8,51.2, ...
       50.9,51.7,52.4,52.9,52.9,52.5,52.0,51.5,51.3,51.4,52.1,53.0, ...
       53.9,54.5,55.0,55.5,56.0,56.4,56.4,56.8,58.0,59.3,59.5,58.7, ...
       57.4,56.1,54.9,54.1,53.5,52.8,51.8,50.8,50.1,49.7,49.8,50.4, ...
       50.9,50.9,50.3,49.7,49.3,49.2,49.3,49.4,49.5,49.6,49.8,50.2];
    NI = 2*24;
    NJ = 10;
    zp_mlat = zeros(NI,NJ);
    k = 1;
    for j=1:NJ
      for i=1:NI
        zp_mlat(i,j) = dat(k);
        k = k + 1;
      end
    end
    ab_mlt = zeros(NI,1);
    for i=1:NI 
      ab_mlt(i)=(i-1)*.5;
    end
  end
  ab_mlat = zeros(NI,1);
        
  if xkp > NJ-1
    xkp=NJ-1;
  end
  kp1=floor(xkp)+1;
  xkp1=floor(xkp)*1.0;
  kp2=kp1+1;
  if kp2 > NJ
    kp2=NJ;
  end

  for i=1:NI 
    ab_mlat(i)=zp_mlat(i,kp1)+(xkp-xkp1)* ...
              (zp_mlat(i,kp2)-zp_mlat(i,kp1));
  end

  cgmlat=-99.99;

  if xmlt < 0.0 || xmlt >= 24.0
    return;
  end

  i1=floor(xmlt/0.5)+1;
  if i1 >= NI
    i1=1;
  end
  i2=i1+1;

  s1=(zp_mlat(i2,kp1)-zp_mlat(i1,kp1))/(ab_mlt(i2)-ab_mlt(i1));
  zmlkp1=zp_mlat(i1,kp1)+(xmlt-ab_mlt(i1))*s1;
  s2=(zp_mlat(i2,kp2)-zp_mlat(i1,kp2))/(ab_mlt(i2)-ab_mlt(i1));
  zmlkp2=zp_mlat(i1,kp2)+(xmlt-ab_mlt(i1))*s2;

  cgmlat=zmlkp1+(xkp-xkp1)*(zmlkp2-zmlkp1);

end

