function [ STORME_AP ] = STORME_AP( context,JDOY,XMLAT,AP )
%STORME_AP EMPIRICAL STORM-E MODEL
%
%      FUNCTION STORME_AP(JDOY,XMLAT,AP)
%CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
%
% EMPIRICAL STORM-E MODEL: COMPUTES A STORM-TO-QUIET RATIO (SQR) FACTOR
% TO ADJUST THE QUIESCENT E-REGION PEAK ELECTRON DENSITY TO ACCOUNT FOR
% ENHANCEMENTS DUE TO GEOMAGNETIC ACTIVITY. THE SQR FACTORS WERE
% COMPUTED FROM NO+ 4.3 UM VOLUME EMISSION RATES DERIVED FROM 
% TIMED/SABER LIMB RADIANCE MEASUREMENTS. THE SABER-DERIVED SQR FACTORS
% WERE FIT TO POWE-LAW IN THE ap INDEX.   
%
% INPUT PARAMETERS:
%
%  JDOY      --- DAY OF YEAR (1-365) 
%  XMLAT     --- MAGNETIC LATITUDE (DEGREES)
%  AP        --- ap INDEX
%
% OUTPUT PARAMETER
% 
%  STORME_AP --- STORM-TO-QUIET RATIO (SQR) TO ADJUST QUIESCENT E-REGION
%                PEAK ELECTRON DENSITY TO ACCOUNT FOR GEOMAGNETIC
%                ENHANCEMENTS. SQR COMPUTED FROM A POWER-LAW FIT
%                IN AP-INDEX: SQR=C1*AP**C2+C3
%
% REFERENCES:
% 
%  (1) Mertens et al. [submitted to JASR, 2011]
%  (2) Fernandez et al. [JASR, Vol. 46, 2010]
%  (3) Mertens et al. [Proc. of SPIE, Vol. 7475, 2009]
%  (4) Mertens et al. [Proc. of SPIE, Vol. 6745, 2007]
%  (5) Mertens et al. [JASR, Vol. 39, 2007] 
% 
% SOFTWARE WRITTEN BY Christopher J. Mertens
%                     NASA Langley Research Center
%                     Atmospheric Sciences Competency
%                     21 Langley Blvd., Mail Stop 401B
%                     Hampton, VA 23681-2199
%
%CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  persistent NMLG NDBD XMLG IDBD C1 C2 C3;
  if isempty(NMLG)
    NMLG=37;
    NDBD=5;
    %  DIMENSION C3(NMLG,NDBD)
    %  DIMENSION XMLG(NMLG),IDBD(NDBD),C1(NMLG,NDBD),C2(NMLG,NDBD)

    %  COMMON /iounit/konsol

    XMLG = [-90.0,-85.0,-80.0,-75.0,-70.0,-65.0,-60.0,-55.0,-50.0, ...
                -45.0,-40.0,-35.0,-30.0,-25.0,-20.0,-15.0,-10.0,-5.0, ...
                0.0,5.0,10.0,15.0,20.0,25.0,30.0,35.0,40.0,45.0,50.0, ...
                55.0,60.0,65.0,70.0,75.0,80.0,85.0,90.0];
    IDBD = [79,171,264,354,366];
    %
    % JDOY        | 0-79 ||80-171||172-264||265-354||355-365|
    %
    C1 = ...
           [ 0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ... -90.0 
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ... -85.0
             0.00000, 0.00508, 0.17360, 0.00000, 0.00000; ... -80.0
             0.00000, 0.31576, 0.31498, 0.00000, 0.00000; ... -75.0
             0.00000, 0.39217, 0.40121, 0.00000, 0.00000; ... -70.0
             0.00000, 0.32634, 0.30179, 0.00000, 0.00000; ... -65.0
             0.11573, 0.06211, 0.16230, 0.20233, 0.11573; ... -60.0
             0.00526, 0.00013, 0.00204, 0.04965, 0.00526; ... -55.0
             0.00011, 0.00013, 0.00018, 0.00040, 0.00011; ... -50.0
             0.00001, 0.00002, 0.00040, 0.00001, 0.00001; ... -45.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ... -40.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ... -35.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ... -30.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ... -25.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ... -20.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ... -15.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ... -10.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ...  -5.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ...   0.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ...   5.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ...  10.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ...  15.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ...  20.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ...  25.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ...  30.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ...  35.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ...  40.0
             0.00000,-0.02738, 0.00004, 0.00001, 0.00000; ...  45.0
             0.00001,-0.00022, 0.00001, 0.00004, 0.00001; ...  50.0
             0.00126, 0.00062, 0.00011, 0.00345, 0.00126; ...  55.0
             0.14923, 0.05483, 0.07113, 0.18282, 0.14923; ...  60.0
             0.37361, 0.00000, 0.00000, 0.44592, 0.37361; ...  65.0
             0.27792, 0.00000, 0.00000, 0.00804, 0.27792; ...  70.0
             0.06445, 0.00000, 0.00000, 0.10315, 0.06445; ...  75.0
             0.00149, 0.00000, 0.00000, 0.00073, 0.00149; ...  80.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ...  85.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000]; %  90.0  
    %    ((C1(IML,IDB),IDB=1,NDBD),IML=1,NMLG)

    %
    % JDOY        | 0-79 ||80-171||172-264||265-354||355-365|
    %
    C2 = ...
           [ 0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ... -90.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ... -85.0
             0.00000, 1.00000, 0.32415, 0.00000, 0.00000; ... -80.0
             0.00000, 0.36538, 0.35455, 0.00000, 0.00000; ... -75.0
             0.00000, 0.41287, 0.38062, 0.00000, 0.00000; ... -70.0
             0.00000, 0.52224, 0.52810, 0.00000, 0.00000; ... -65.0
             0.73025, 0.90723, 0.68107, 0.64815, 0.73025; ... -60.0
             1.29410, 2.06038, 1.47332, 0.84843, 1.29410; ... -55.0
             1.79442, 1.77511, 1.59906, 1.59141, 1.79442; ... -50.0
             1.84434, 1.70607, 1.03056, 1.92168, 1.84434; ... -45.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ... -40.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ... -35.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ... -30.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ... -25.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ... -20.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ... -15.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ... -10.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ...  -5.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ...   0.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ...   5.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ...  10.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ...  15.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ...  20.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ...  25.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ...  30.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ...  35.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ...  40.0
             2.04797, 0.27034, 1.00000, 1.90792, 2.04797; ...  45.0
             2.24520, 1.00000, 1.88721, 2.01535, 2.24520; ...  50.0
             1.56493, 1.52468, 1.93389, 1.38532, 1.56493; ...  55.0
             0.71442, 0.87492, 0.78890, 0.66828, 0.71442; ...  60.0
             0.53546, 0.00000, 0.00000, 0.42597, 0.53546; ...  65.0
             0.48647, 0.00000, 0.00000, 1.00000, 0.48647; ...  70.0
             0.67340, 0.00000, 0.00000, 0.36809, 0.67340; ...  75.0
             1.44025, 0.00000, 0.00000, 1.13529, 1.44025; ...  80.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000; ...  85.0
             0.00000, 0.00000, 0.00000, 0.00000, 0.00000]; %  90.0 
    %
    % JDOY        | 0-79 ||80-171||172-264||265-354||355-365|
    %
    C3 = ...
           [ 1.00000, 1.00000, 1.00000, 1.00000, 1.00000; ... -90.0
             1.00000, 1.00000, 1.00000, 1.00000, 1.00000; ... -85.0
             1.00000, 1.03132, 0.76703, 1.00000, 1.00000; ... -80.0
             1.00000, 0.57588, 0.56324, 1.00000, 1.00000; ... -75.0
             1.00000, 0.41370, 0.38549, 1.00000, 1.00000; ... -70.0
             1.00000, 0.51704, 0.50217, 1.00000, 1.00000; ... -65.0
             0.55236, 0.80162, 0.60824, 0.46999, 0.55236; ... -60.0
             0.90923, 0.99688, 0.96752, 0.67312, 0.90923; ... -55.0
             0.99338, 0.98486, 0.99503, 0.87473, 0.99338; ... -50.0
             1.00031, 1.00369, 1.00225, 0.91242, 1.00031; ... -45.0
             1.00000, 1.00000, 1.00000, 1.00000, 1.00000; ... -40.0
             1.00000, 1.00000, 1.00000, 1.00000, 1.00000; ... -35.0
             1.00000, 1.00000, 1.00000, 1.00000, 1.00000; ... -30.0
             1.00000, 1.00000, 1.00000, 1.00000, 1.00000; ... -25.0
             1.00000, 1.00000, 1.00000, 1.00000, 1.00000; ... -20.0
             1.00000, 1.00000, 1.00000, 1.00000, 1.00000; ... -15.0
             1.00000, 1.00000, 1.00000, 1.00000, 1.00000; ... -10.0
             1.00000, 1.00000, 1.00000, 1.00000, 1.00000; ...  -5.0
             1.00000, 1.00000, 1.00000, 1.00000, 1.00000; ...   0.0
             1.00000, 1.00000, 1.00000, 1.00000, 1.00000; ...   5.0
             1.00000, 1.00000, 1.00000, 1.00000, 1.00000; ...  10.0
             1.00000, 1.00000, 1.00000, 1.00000, 1.00000; ...  15.0
             1.00000, 1.00000, 1.00000, 1.00000, 1.00000; ...  20.0
             1.00000, 1.00000, 1.00000, 1.00000, 1.00000; ...  25.0
             1.00000, 1.00000, 1.00000, 1.00000, 1.00000; ...  30.0
             1.00000, 1.00000, 1.00000, 1.00000, 1.00000; ...  35.0
             1.00000, 1.00000, 1.00000, 1.00000, 1.00000; ...  40.0
             1.04797, 1.02319, 0.98581, 1.01457, 1.04797; ...  45.0
             1.03332, 0.97016, 0.97807, 0.99044, 1.03332; ...  50.0
             1.00633, 0.94822, 0.96340, 0.95363, 1.00633; ...  55.0
             0.67902, 0.71540, 0.70230, 0.60821, 0.67902; ...  60.0
             0.35017, 1.00000, 1.00000, 0.51033, 0.35017; ...  65.0
             0.63358, 1.00000, 1.00000, 1.37782, 0.63358; ...  70.0
             0.85724, 1.00000, 1.00000, 0.91942, 0.85724; ...  75.0
             0.92703, 1.00000, 1.00000, 1.00502, 0.92703; ...  80.0
             1.00000, 1.00000, 1.00000, 1.00000, 1.00000; ...  85.0
             1.00000, 1.00000, 1.00000, 1.00000, 1.00000]; %  90.0
  end
  %
  % ... Find Season-Averaged Coefficient Index 
  %
  IDXS=0;
  if JDOY <= IDBD(1)
    IDXS=1;
  end
  for IS=2:NDBD
    if (JDOY > IDBD(IS-1)) && (JDOY <= IDBD(IS))
      IDXS=IS;
    end
  end
  if IDXS == 0
    if context.KONSOL > 0
      fprintf(context.KONSOL,'ERROR IN STORME_AP: PROBLEM FINDING SEASON-AVERAGED COEFFICIENT\n');
      fprintf(context.KONSOL,'   INPUT DAY OF YEAR = %g\n',JDOY);
    end
    STORME_AP=-5.0;
  else
    %
    % ... Find Magnetic Latitude Coefficient Index
    % 
    IDXL=0;
    DELG=abs(XMLG(1)-XMLG(2));
    DELD=DELG/2.0;
    YMP=XMLG(1)+DELD;
    YMM=XMLG(NMLG)-DELD;
    if (XMLAT >= XMLG(1)) && (XMLAT <= YMP)
      IDXL=1;
    end
    if (XMLAT > YMM) && (XMLAT <= XMLG(NMLG))
      IDXL=NMLG;
    end
    for IL=2:NMLG-1
      YMP=XMLG(IL)+DELD;
      YMM=XMLG(IL)-DELD;
      if (XMLAT > YMM) && (XMLAT <= YMP)
        IDXL=IL;
      end
    end
    if IDXL == 0
      if context.KONSOL > 0
        fprintf(context.KONSOL,'ERROR IN STORME_AP: PROBLEM FINDING MAGNETIC LATITUDE COEFFICIENT\n');
        fprintf(context.KONSOL,'     INPUT MAGNETIC LATITUDE(DEGREES) = %g\n',XMLAT);
      end
      STORME_AP=-5.0;
    else
      %
      % ... COMPUTE E-REGION ELECTRON DENSITY GEOMAGNETIC STORM ENHANCEMET
      % ... FACTOR (i.e., THE STORM-TO-QUIET RATIO (SQR)) 
      %
      SQR=C1(IDXL,IDXS)*AP^(C2(IDXL,IDXS))+C3(IDXL,IDXS);
      if SQR < 1.0
        SQR=1.0;
      end
      STORME_AP=SQR;
    end
  end

end

