function [ B,context ] = SHAB1D( context, FLAT,FLON,T,RZ )
%SHAB1D COMPUTES THE HOURLY VALUES OF B1 FROM A SET OF SH COEFFICIENTS
%
%      SUBROUTINE  SHAB1D (FLAT,FLON,T,RZ,B)
%-------------------------------------------------------------------
%	COMPUTES THE HOURLY VALUES OF B1 FROM A SET OF SH COEFFICIENTS
%	IN A POINT OF A GIVEN GEOCENTRIC LATITUDE AND LONGITUDE
%	OF THE EARTH'S SURFACE FOR A GIVEN MONTH AND A GIVEN SUSPOT NUMER
%
%   PARAMETERS ARE THE SAME AS IN SHAMDB0D, EXCEPT:
%		FLAT	Geographic latitude
%		B		=B1
%
%	***** PARAMS & COEFFS IN DATA SENTENCES *****
%-------------------------------------------------------------------

%      PARAMETER   (IBO=0,JBO=1,KDIM=6,LDIM=4,L=-1)
%      DIMENSION   FN(0:KDIM,0:KDIM), CONST(0:KDIM,0:KDIM)
%      DIMENSION   GNM(0:KDIM,0:KDIM,1-IBO-JBO:LDIM),
%     *			HNM(0:KDIM,0:KDIM,1-IBO-JBO:LDIM)
%	  DIMENSION   GANM(0:KDIM,0:KDIM,1-IBO-JBO:LDIM),
%     *			HANM(0:KDIM,0:KDIM,1-IBO-JBO:LDIM),
%     *            GBNM(0:KDIM,0:KDIM,1-IBO-JBO:LDIM),
%     *			HBNM(0:KDIM,0:KDIM,1-IBO-JBO:LDIM)
%      DIMENSION   BINT(0:KDIM,0:KDIM,1-IBO-JBO:LDIM),
%     *            BEXT(0:KDIM,0:KDIM,1-IBO-JBO:LDIM)
%      CHARACTER*1 IE       
%      COMMON  BINT,BEXT,RE,TZERO,IFIT,IB,KINT,LINT,KEXT,
%     *              LEXT,KMAX,FN
%c      ,CONST
  persistent IBO JBO KDIM LDIM L RE TZERO IFIT IB KINT ...
    LINT KEXT LEXT IE GANM GBNM HANM HBNM KMAX KT;
  if isempty(IBO)
    IBO=0;
    JBO=1;
    KDIM=6;
    LDIM=4;
    L=-1;
	  %ALT = 300.;

    %THETA = 180.;
    %RE = 6371.2;
    RE=WGS84.equatorialRadius/1000.0; % km
    TZERO = 1.0;
    IFIT = -1;
    %ICEN = 0;
    %IREF = 0;
    IB = 2;
    KINT = 6;
    LINT = 4;
    KEXT = 0;
    LEXT = -1;

%   	cnst = ...
%       	 [1.,1.,1.,1.,1.73205,0.866025,1.,2.44949,1.93649,0.7905691,1., ...
%       	  3.16228,3.35410,2.09165,0.739510,1.,3.87298,5.12348, ...
%       	  4.18330,2.21853,0.701561,1.,4.58258,7.24569,7.24569,4.96078, ...
%             2.32681,0.671693];
%     CONST = zeros(KDIM+1,KDIM+1);
%     K = 1;
%     for N=0:KDIM
%       for M=0:N
%         CONST(N+1,M+1) = cnst(K);
%         %CONST(M+1,N+1) = cnst(K);
%         K = K + 1;
%       end
%     end

	  IE = 'I';

    dat = ...
         [1.156, 0.039,  0.000, 0.000, ...
          1.725,-0.053,  0.000, 0.000,  1.097,-0.032,  0.000, 0.000, ...
         -0.579, 0.019,  0.000, 0.000,  0.265,-0.010,  0.000, 0.000, ...
         -2.895, 0.023,  0.000, 0.000, ...
          3.269,-0.025,  0.000, 0.000,  2.278,-0.015,  0.000, 0.000, ...
         -1.789, 0.009,  0.000, 0.000,  0.653,-0.005,  0.000, 0.000, ...
         -3.240, 0.052, -2.645, 0.030, ...
          4.404,-0.062,  3.283,-0.038,  2.827,-0.038,  2.181,-0.019, ...
         -1.496, 0.020, -1.143, 0.011,  0.688,-0.011,  0.512,-0.008, ...
         -0.023,-0.025,  0.000, 0.000, ...
         -0.370, 0.031,  0.000, 0.000, -0.385, 0.017,  0.000, 0.000, ...
          0.150,-0.009,  0.000, 0.000,  0.019, 0.007,  0.000, 0.000, ...
          1.704, 0.006, -1.766, 0.022, ...
         -2.115,-0.007,  2.309,-0.028, -1.610,-0.003,  1.314,-0.015, ...
          0.806, 0.002, -0.873, 0.010, -0.249,-0.002,  0.445,-0.005, ...
          0.256,-0.009, -2.959, 0.023, ...
         -0.662, 0.013,  3.643,-0.030, -0.208, 0.001,  2.208,-0.015, ...
          0.144,-0.001, -1.326, 0.010, -0.205, 0.005,  0.678,-0.007, ...
          1.730,-0.030,  0.000, 0.000, ...
         -2.072, 0.035,  0.000, 0.000, -1.027, 0.016,  0.000, 0.000, ...
          0.946,-0.008,  0.000, 0.000, -0.644, 0.009,  0.000, 0.000, ...
          3.925,-0.060,  1.613,-0.015, ...
         -4.790, 0.070, -2.021, 0.019, -3.293, 0.048, -1.273, 0.006, ...
          1.829,-0.026,  0.797,-0.005, -0.767, 0.011, -0.395, 0.006, ...
         -1.988, 0.027, -0.761, 0.003, ...
          2.258,-0.031,  0.978,-0.004,  1.772,-0.026,  0.459,-0.003, ...
         -0.813, 0.014, -0.257, 0.004,  0.173,-0.003,  0.223,-0.001, ...
          0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000, ...
          0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000, ...
         -1.356, 0.011,  0.000, 0.000, ...
          0.079,-0.003,  0.000, 0.000,  0.415,-0.012,  0.000, 0.000, ...
         -0.249, 0.004,  0.000, 0.000, -0.087, 0.004,  0.000, 0.000, ...
         -1.155,-0.012,  0.261,-0.026, ...
          1.619, 0.011, -0.217, 0.032,  0.989, 0.009, -0.361, 0.019, ...
         -0.745,-0.001,  0.009,-0.009,  0.347, 0.000,  0.178, 0.005, ...
          4.672,-0.032,  0.562, 0.017, ...
         -5.868, 0.040, -0.850,-0.020, -3.798, 0.026, -0.145,-0.020, ...
          2.079,-0.014,  0.160, 0.008, -0.943, 0.007, -0.417, 0.001, ...
          0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000, ...
          0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000, ...
          0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000, ...
          0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000, ...
          2.477, 0.000,  0.000, 0.000, ...
         -1.815,-0.011,  0.000, 0.000, -1.571, 0.002,  0.000, 0.000, ...
          0.551, 0.004,  0.000, 0.000, -0.044,-0.008,  0.000, 0.000, ...
          2.160,-0.010, -0.305, 0.012, ...
         -2.618, 0.010,  0.529,-0.016, -1.782, 0.006,  0.634,-0.007, ...
          0.976,-0.003, -0.449, 0.006, -0.331, 0.001, -0.004,-0.004, ...
         -0.394, 0.002,  0.851,-0.020, ...
          0.359,-0.003, -1.051, 0.024,  0.357, 0.002, -0.239, 0.012, ...
          0.005,-0.001,  0.210,-0.008, -0.028,-0.003, -0.322, 0.005, ...
          0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000, ...
          0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000, ...
          0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000, ...
          0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000, ...
          0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000, ...
          0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000, ...
         -6.760, 0.064,  0.000, 0.000, ...
          7.700,-0.073,  0.000, 0.000,  5.394,-0.054,  0.000, 0.000, ...
         -2.788, 0.026,  0.000, 0.000,  0.923,-0.007,  0.000, 0.000, ...
         -2.328, 0.024, -0.463, 0.020, ...
          2.923,-0.027,  0.490,-0.025,  1.768,-0.019,  0.711,-0.017, ...
         -1.068, 0.009, -0.363, 0.010,  0.596,-0.004, -0.073,-0.004, ...
         -1.911, 0.016, -4.519, 0.041, ...
          2.644,-0.024,  5.569,-0.050,  1.287,-0.009,  3.707,-0.031, ...
         -0.894, 0.007, -2.121, 0.019,  0.669,-0.007,  0.933,-0.010, ...
          0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000, ...
          0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000, ...
          0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000, ...
          0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000, ...
          0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000, ...
          0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000, ...
          0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000, ...
          0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000,0.000];
    GANM = zeros(KDIM+1,KDIM+1,LDIM-(1-IBO-JBO)+1);
    GBNM = zeros(KDIM+1,KDIM+1,LDIM-(1-IBO-JBO)+1);
    HANM = zeros(KDIM+1,KDIM+1,LDIM-(1-IBO-JBO)+1);
    HBNM = zeros(KDIM+1,KDIM+1,LDIM-(1-IBO-JBO)+1);
    K = 1;
    for N=0:KDIM
      for M=0:N
        for J=1-IBO-JBO:LDIM
          GANM(N+1,M+1,J-(1-IBO-JBO)+1) = dat(K);
          %GANM(M+1,N+1,J-(1-IBO-JBO)+1) = dat(K);
          GBNM(N+1,M+1,J-(1-IBO-JBO)+1) = dat(K+1);
          %GBNM(M+1,N+1,J-(1-IBO-JBO)+1) = dat(K+1);
          HANM(N+1,M+1,J-(1-IBO-JBO)+1) = dat(K+2);
          %HANM(M+1,N+1,J-(1-IBO-JBO)+1) = dat(K+2);
          HBNM(N+1,M+1,J-(1-IBO-JBO)+1) = dat(K+3);
          %HBNM(M+1,N+1,J-(1-IBO-JBO)+1) = dat(K+3);
          K = K + 4;
        end
      end
    end

    KMAX = max(KINT,KEXT);
    if KMAX > KDIM
      % stop
      fprintf(1,'ERROR - bad KINT(%d) or KEXT(%d) in SHAB1D (%d > %d)\n',KINT,KEXT,KMAX,KDIM);
      return;
    end
    KT = max(LINT,LEXT);
    if KT > LDIM
      % stop
      fprintf(1,'ERROR - bad LINT(%d) or LEXT(%d) in SHAB1D (%d > %d)\n',LINT,LEXT,KT,LDIM);
      return;
    end
  end
  context.RE = RE;
  context.TZERO = TZERO;
  context.IFIT = IFIT;
  context.IB = IB;
  context.KINT = KINT;
  context.KINT = KINT;
  context.LINT = LINT;
  context.KEXT = KEXT;
  context.LEXT = LEXT;
  context.KMAX = KMAX;
  context.FN = zeros(KDIM+1,KDIM+1);
  context.BINT = zeros(KDIM+1,KDIM+1,LDIM-(1-IBO-JBO)+1);
  context.BEXT = zeros(KDIM+1,KDIM+1,LDIM-(1-IBO-JBO)+1);
  GNM = zeros(KDIM+1,KDIM+1,LDIM-(1-IBO-JBO)+1);
  HNM = zeros(KDIM+1,KDIM+1,LDIM-(1-IBO-JBO)+1);

	for N=0:KMAX
    for M=0:N

      for J=1-IBO-JBO:KT
        GNM(N+1,M+1,J-(1-IBO-JBO)+1)=GANM(N+1,M+1,J-(1-IBO-JBO)+1) ...
                                    +GBNM(N+1,M+1,J-(1-IBO-JBO)+1)*RZ;
        HNM(N+1,M+1,J-(1-IBO-JBO)+1)=HANM(N+1,M+1,J-(1-IBO-JBO)+1) ...
                                    +HBNM(N+1,M+1,J-(1-IBO-JBO)+1)*RZ;
      end
 
      if IE == 'I'
        if N > KINT
          continue;
        end
        LJ = LINT;
      else
        if N > KEXT
          continue;
        end
        LJ = LEXT;
      end

      context.FN(N+1,M+1) = N;

      if M <= 0
        for J=1-IBO-JBO:KT
          if IE == 'I'
            context.BINT(N+1,M+1,J-(1-IBO-JBO)+1)   = GNM(N+1,M+1,J-(1-IBO-JBO)+1);
          else
            context.BEXT(N+1,M+1,J-(1-IBO-JBO)+1)   = GNM(N+1,M+1,J-(1-IBO-JBO)+1);
          end
        end
      else
        for J=1-IBO-JBO:LJ
          if IE == 'I'
            context.BINT(N+1,M+1,J-(1-IBO-JBO)+1)   = GNM(N+1,M+1,J-(1-IBO-JBO)+1);
            context.BINT(M-1+1,N+1,J-(1-IBO-JBO)+1) = HNM(N+1,M+1,J-(1-IBO-JBO)+1);
          else
            context.BEXT(N+1,M+1,J-(1-IBO-JBO)+1)   = GNM(N+1,M+1,J-(1-IBO-JBO)+1);
            context.BEXT(M-1+1,N+1,J-(1-IBO-JBO)+1) = HNM(N+1,M+1,J-(1-IBO-JBO)+1);
          end
        end
      end
    end
  end
  %     **********************************************************
  %     SYNTHESIZES THE VALUE OF B1 FROM THE MODEL	
  %     **********************************************************
  [~,~,B,context] = SCHNEVPD(context,RZ,FLAT,FLON,0,T,L);

end

