function [BX,BY,BZ] = T01_01(IOPT,PARMOD,PS,X,Y,Z);
% function [BX,BY,BZ] = T01(IOPT,PARMOD,PS,X,Y,Z);
% Tsyganenko's External Field Model, 2001 Version (T01)
% Translated from original FORTRAN March 27, 2003
% By Paul O'Brien (original by N.A. Tsyganenko)
% Paul.OBrien@aero.org (Nikolai.Tsyganenko@gsfc.nasa.gov)
%
% All subroutines enclosed here as subfunctions
%
% Updates:
%   10 September, 2004 - D. Berube found some initialization errors in EXTALL and DIPOLE
%   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%      SUBROUTINE T01_01 (IOPT,PARMOD,PS,X,Y,Z,BX,BY,BZ)
% c
% c     RELEASE DATE OF THIS VERSION:   AUGUST 8, 2001.
% C
% c--------------------------------------------------------------------
% C   A DATA-BASED MODEL OF THE EXTERNAL (I.E., WITHOUT EARTH'S CONTRIBUTION) PART OF THE
% C   MAGNETOSPHERIC MAGNETIC FIELD, CALIBRATED BY
% C    (1) SOLAR WIND PRESSURE PDYN (NANOPASCALS),
% C    (2) DST (NANOTESLA),
% C    (3) BYIMF,
% C    (4) BZIMF (NANOTESLA)
% C    (5) G1-INDEX
% C    (6) G2-INDEX  (SEE TSYGANENKO [2001] FOR AN EXACT DEFINITION OF THESE TWO INDICES)
% 
% c   THESE INPUT PARAMETERS SHOULD BE PLACED IN THE FIRST 6 ELEMENTS
% c   OF THE ARRAY PARMOD(10).
% C
% C   THE REST OF THE INPUT VARIABLES ARE: THE GEODIPOLE TILT ANGLE PS (RADIANS),
% C     AND   X,Y,Z -  GSM POSITION (RE)
% C
% c   IOPT  IS JUST A DUMMY INPUT PARAMETER, NECESSARY TO MAKE THIS SUBROUTINE
% C   COMPATIBLE WITH THE TRACING SOFTWARE PACKAGE (GEOPACK). IN THIS MODEL
% C   IT DOES NOT AFFECT THE OUTPUT FIELD.
% c
% C*******************************************************************************************
% c** ATTENTION:  THE MODEL IS BASED ON DATA TAKEN SUNWARD FROM X=-15Re, AND HENCE BECOMES   *
% C**              INVALID AT LARGER TAILWARD DISTANCES !!!                                  *
% C*******************************************************************************************
% C
% c   OUTPUT:  GSM COMPONENTS OF THE EXTERNAL MAGNETIC FIELD (BX,BY,BZ, nanotesla)
% C            COMPUTED AS A SUM OF CONTRIBUTIONS FROM PRINCIPAL FIELD SOURCES
% C
% c  (C) Copr. 2001, Nikolai A. Tsyganenko, USRA, Code 690.2, NASA GSFC
% c      Greenbelt, MD 20771, USA
% c
% C                            REFERENCE:
% C
% C    N. A. Tsyganenko, A new data-based model of the near magnetosphere magnetic field:
% c       1. Mathematical structure.
% c       2. Parameterization and fitting to observations.
% c
% c             (submitted to JGR, July 2001; available online in the PDF format
% c              from anonymous ftp-area www-istp.gsfc.nasa.gov,  /pub/kolya/T01)
% c
% c----------------------------------------------------------------------
% c
% c
% c This is a sample main program, calculating the T01 model field produced by the
% c external (magnetospheric) sources of the geomagnetic field at a specified point
% c of space {XGSM,YGSM,ZGSM}.  The earth's dipole tilt angle PS should be specified
% c in radians, the solar wind ram pressure PDYN in nanoPascals, Dst index, IMF By
% c and Bz components in nT. The two IMF-related indices G1 and G2 take into account
% c the IMF and solar wind conditions during the preceding 1-hour interval; their
% c exact definition is given in the paper "A new data-based model of the near
% c magnetosphere magnetic field. 2. Parameterization and fitting to observations".
% c The paper is available online from anonymous ftp-area www-istp.gsfc.nasa.gov,
% c /pub/kolya/T01.
% c
% % function test01
% % parmod = repmat(0,10,1);
% % parmod(1)=2;   %PDYN
% % parmod(2)=-20; %DST
% % parmod(3)=-1;  %BYIMF
% % parmod(4)=-2;  %BZIMF
% % parmod(5)=10;   %G1
% % parmod(6)=10;   %G2
% % iopt = 3;
% % ps = 0.5;
% % x = 1.5;
% % y = 2.5;
% % z = 3.5;
% % [bx,by,bz] = T01_01(iopt,parmod,ps,x,y,z);
% % disp([bx,by,bz]);



%      REAL PARMOD(10),PS,X,Y,Z,BX,BY,BZ
%      REAL*8 A(43),PDYN,DST_AST,BYIMF,BZIMF,G1,G2,PSS,XX,YY,ZZ,
%     *  BXCF,BYCF,BZCF,BXT1,BYT1,BZT1,BXT2,BYT2,BZT2,
%     *  BXSRC,BYSRC,BZSRC,BXPRC,BYPRC,BZPRC, BXR11,BYR11,BZR11,
%     *  BXR12,BYR12,BZR12,BXR21,BYR21,BZR21,BXR22,BYR22,BZR22,HXIMF,
%     *  HYIMF,HZIMF,BBX,BBY,BBZ
% C
A = [1.00000,2.48341,.58315,.31917,-.08796,-1.17266,3.57478, ...
        -.06143,-.01113,.70924,-.01675,-.46056,-.87754,-.03025,.18933, ...
        .28089,.16636,-.02932,.02592,-.23537,-.07659,.09117,-.02492, ...
        .06816,.55417,.68918,-.04604,2.33521,3.90147,1.28978,.03139, ...
        .98751,.21824,41.60182,1.12761,.01376,1.02751,.02969,.15790, ...
        8.94335,28.31280,1.24364,.38013];
% C
if (X<-20.),
    disp(sprintf(['  ATTENTION:  THE MODEL IS VALID SUNWARD FROM X=-15 Re ONLY,\n',...
            '              WHILE YOU ARE TRYING TO USE IT AT X=%g'], X));
%    pause;
end
% C
PDYN=PARMOD(1);
DST_AST=PARMOD(2)*0.8-13.*sqrt(PDYN);
BYIMF=PARMOD(3);
BZIMF=PARMOD(4);
% C
G1=PARMOD(5);
G2=PARMOD(6);
PSS=PS;
XX=X;
YY=Y;
ZZ=Z;
% C

[BXCF,BYCF,BZCF,BXT1,BYT1,BZT1,BXT2,BYT2,BZT2, ...
        BXSRC,BYSRC,BZSRC,BXPRC,BYPRC,BZPRC, BXR11,BYR11,BZR11, ...
        BXR12,BYR12,BZR12,BXR21,BYR21,BZR21,BXR22,BYR22,BZR22,HXIMF,...
        HYIMF,HZIMF,BBX,BBY,BBZ] = T01_EXTALL(0,0,0,0,A,43,PDYN,DST_AST, ...
    BYIMF,BZIMF,G1,G2,PSS,XX,YY,ZZ);
% C
BX=BBX;
BY=BBY;
BZ=BBZ;
% C
% end of function T01_01

% C
% C================================================================
function [BXCF,BYCF,BZCF,BXT1,BYT1,BZT1,BXT2,BYT2,BZT2, ...
        BXSRC,BYSRC,BZSRC,BXPRC,BYPRC,BZPRC, BXR11,BYR11,BZR11, ...
        BXR12,BYR12,BZR12,BXR21,BYR21,BZR21,BXR22,BYR22,BZR22,HXIMF, ...
        HYIMF,HZIMF,BX,BY,BZ] = T01_EXTALL(IOPGEN,IOPT,IOPB,IOPR,A,NTOT, ...
    PDYN,DST,BYIMF,BZIMF,VBIMF1,VBIMF2,PS,X,Y,Z);
%       SUBROUTINE EXTALL (IOPGEN,IOPT,IOPB,IOPR,A,NTOT,
%      *  PDYN,DST,BYIMF,BZIMF,VBIMF1,VBIMF2,PS,X,Y,Z,
%      *  BXCF,BYCF,BZCF,BXT1,BYT1,BZT1,BXT2,BYT2,BZT2,
%      *  BXSRC,BYSRC,BZSRC,BXPRC,BYPRC,BZPRC, BXR11,BYR11,BZR11,
%      *  BXR12,BYR12,BZR12,BXR21,BYR21,BZR21,BXR22,BYR22,BZR22,HXIMF,
%      *  HYIMF,HZIMF,BX,BY,BZ)
% C
% C   IOPGEN - GENERAL OPTION FLAG:  IOPGEN=0 - CALCULATE TOTAL FIELD
% C                                  IOPGEN=1 - DIPOLE SHIELDING ONLY
% C                                  IOPGEN=2 - TAIL FIELD ONLY
% C                                  IOPGEN=3 - BIRKELAND FIELD ONLY
% C                                  IOPGEN=4 - RING CURRENT FIELD ONLY
% C                                  IOPGEN=5 - INTERCONNECTION FIELD ONLY
% C
% C   IOPT -  TAIL FIELD FLAG:       IOPT=0  -  BOTH MODES
% C                                  IOPT=1  -  MODE 1 ONLY
% C                                  IOPT=2  -  MODE 2 ONLY
% C
% C   IOPB -  BIRKELAND FIELD FLAG:  IOPB=0  -  ALL 4 TERMS
% C                                  IOPB=1  -  REGION 1, MODES 1 AND 2
% C                                  IOPB=2  -  REGION 2, MODES 1 AND 2
% C
% C   IOPR -  RING CURRENT FLAG:     IOPR=0  -  BOTH SRC AND PRC
% C                                  IOPR=1  -  SRC ONLY
% C                                  IOPR=2  -  PRC ONLY
% C
%      IMPLICIT  REAL * 8  (A - H, O - Z)
% C
%      DIMENSION A(NTOT)
% C
%       COMMON /TAIL/ DXSHIFT1,DXSHIFT2,D,DELTADY  %! THE COMMON BLOCKS FORWARD NONLINEAR PARAMETERS
%       COMMON /BIRKPAR/ XKAPPA1,XKAPPA2
%       COMMON /RCPAR/ SC_SY,SC_AS,PHI
%       COMMON /G/ G
%       COMMON /RH0/ RH0

global T01; % common blocks G and RHO are scalars, not structures

[BXCF,BYCF,BZCF,BXT1,BYT1,BZT1,BXT2,BYT2,BZT2, ...
        BXSRC,BYSRC,BZSRC,BXPRC,BYPRC,BZPRC, BXR11,BYR11,BZR11, ...
        BXR12,BYR12,BZR12,BXR21,BYR21,BZR21,BXR22,BYR22,BZR22,HXIMF, ...
        HYIMF,HZIMF,BX,BY,BZ]  = deal(nan); % initialize all outputs to nan
% C
%      DATA A0_A,A0_S0,A0_X0 /34.586D0,1.1960D0,3.4397D0/   %!   SHUE ET AL. PARAMETERS
A0_A = 34.586D0;
A0_S0 = 1.1960D0;
A0_X0 = 3.4397D0;

DSIG = 0.003D0;
if isempty(T01) | ~isfield(T01,'RH0'),
    T01.RH0 = 8.0D0;
end
RH2 = -5.2D0;
% c
XAPPA=(PDYN/2.)^A(39);   %!  NOW THIS IS A VARIABLE PARAMETER
T01.RH0=A(40);
T01.G=A(41);

XAPPA3=XAPPA^3;

XX=X*XAPPA;
YY=Y*XAPPA;
ZZ=Z*XAPPA;
% C
SPS=sin(PS);
% c
X0=A0_X0/XAPPA;
AM=A0_A/XAPPA;
S0=A0_S0;
% c
BPERP=sqrt(BYIMF^2+BZIMF^2);
% C
% C   CALCULATE THE IMF CLOCK ANGLE:
% C
if (BYIMF==0.D0) & (BZIMF==0.D0),
    THETA=0.D0;
else
    THETA=atan2(BYIMF,BZIMF);
    if (THETA<=0.D0), THETA=THETA+6.283185307D0; end
end
% c
CT=cos(THETA);
ST=sin(THETA);
YS=Y*CT-Z*ST;
ZS=Z*CT+Y*ST;

STHETAH=sin(THETA/2.)^2;
% C
% C  CALCULATE "IMF" COMPONENTS OUTSIDE THE MAGNETOPAUSE LAYER (HENCE BEGIN WITH "O")
% C  THEY ARE NEEDED ONLY IF THE POINT (X,Y,Z) IS WITHIN THE TRANSITION MAGNETOPAUSE LAYER
% C  OR OUTSIDE THE MAGNETOSPHERE:
% C
FACTIMF=A(24)+A(25)*STHETAH;
% c
OIMFX=0.D0;
OIMFY=BYIMF*FACTIMF;
OIMFZ=BZIMF*FACTIMF;
% c
R=sqrt(X^2+Y^2+Z^2);
XSS=X;
ZSS=Z;

while 1,
    XSOLD=XSS;      %!   BEGIN ITERATIVE SEARCH OF UNWARPED COORDS (TO FIND SIGMA)
    ZSOLD=ZSS;
    
    RH=T01.RH0+RH2*(ZSS/R)^2;
    SINPSAS=SPS/(1.D0+(R/RH)^3)^0.33333333D0;
    COSPSAS=sqrt(1.D0-SINPSAS^2);
    ZSS=X*SINPSAS+Z*COSPSAS;
    XSS=X*COSPSAS-Z*SINPSAS;
    DD=abs(XSS-XSOLD)+abs(ZSS-ZSOLD);
    if (DD<=1.D-6), break; end
end
% C                                END OF ITERATIVE SEARCH
RHO2=Y^2+ZSS^2;
ASQ=AM^2;
XMXM=AM+XSS-X0;
if (XMXM<0.), XMXM=0; end; %! THE BOUNDARY IS A CYLINDER TAILWARD OF X=X0-AM
AXX0=XMXM^2;
ARO=ASQ+RHO2;
SIGMA=sqrt((ARO+AXX0+sqrt((ARO+AXX0)^2-4.*ASQ*AXX0))/(2.*ASQ));
% C
% C   NOW, THERE ARE THREE POSSIBLE CASES:
% C    (1) INSIDE THE MAGNETOSPHERE   (SIGMA
% C    (2) IN THE BOUNDARY LAYER
% C    (3) OUTSIDE THE MAGNETOSPHERE AND B.LAYER
% C       FIRST OF ALL, CONSIDER THE CASES (1) AND (2):
% C
% C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if (SIGMA<S0+DSIG),  %!  CASES (1) OR (2); CALCULATE THE MODEL FIELD
    % C                              (WITH THE POTENTIAL "PENETRATED" INTERCONNECTION FIELD):
    % C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    % C
    if (IOPGEN<=1),
        [CFX,CFY,CFZ] = T01_SHLCAR3X3(XX,YY,ZZ,PS);         %!  DIPOLE SHIELDING FIELD
        BXCF=CFX*XAPPA3;
        BYCF=CFY*XAPPA3;
        BZCF=CFZ*XAPPA3;
    else
        BXCF=0.D0;
        BYCF=0.D0;
        BZCF=0.D0;
    end
    
    if (IOPGEN==0) | (IOPGEN==2),
        T01.TAIL.DXSHIFT1=A(26)+A(27)*VBIMF2;
        T01.TAIL.DXSHIFT2=0.D0;
        T01.TAIL.D=A(28);
        T01.TAIL.DELTADY=A(29);
        [BXT1,BYT1,BZT1,BXT2,BYT2,BZT2] = T01_DEFORMED(IOPT,PS,XX,YY,ZZ);  %!  TAIL FIELD (THREE MODES)
    else
        BXT1=0.D0;
        BYT1=0.D0;
        BZT1=0.D0;
        BXT2=0.D0;
        BYT2=0.D0;
        BZT2=0.D0;
    end
    
    if (IOPGEN==0) | (IOPGEN==3),
        T01.BIRKPAR.XKAPPA1=A(35)+A(36)*VBIMF2;
        T01.BIRKPAR.XKAPPA2=A(37)+A(38)*VBIMF2;
        [BXR11,BYR11,BZR11,BXR12,BYR12,BZR12,...
                BXR21,BYR21,BZR21,BXR22,BYR22,BZR22] = ...
            T01_BIRK_TOT(IOPB,PS,XX,YY,ZZ)    ;
        %!   BIRKELAND FIELD (TWO MODES FOR R1 AND TWO MODES FOR R2)
    else
        BXR11=0.D0;
        BYR11=0.D0;
        BZR11=0.D0;
        BXR12=0.D0;
        BYR12=0.D0;
        BZR12=0.D0;
        BXR21=0.D0;
        BYR21=0.D0;
        BZR21=0.D0;
        BXR22=0.D0;
        BYR22=0.D0;
        BZR22=0.D0;
    end
    
    if (IOPGEN==0) | (IOPGEN==4),
        T01.RCPAR.PHI=1.5707963D0*tanh(abs(DST)/A(34));
        ZNAM=abs(DST);
        if (ZNAM<20.D0), ZNAM=20.D0; end
        T01.RCPAR.SC_SY=A(30)*(20.D0/ZNAM)^A(31) *XAPPA;    %!
        T01.RCPAR.SC_AS=A(32)*(20.D0/ZNAM)^A(33) *XAPPA;
        [BXSRC,BYSRC,BZSRC,BXPRC,BYPRC,BZPRC] = ...
            T01_FULL_RC(IOPR,PS,XX,YY,ZZ);  %!  SHIELDED RING CURRENT (SRC AND PRC)
    else
        BXSRC=0.D0;
        BYSRC=0.D0;
        BZSRC=0.D0;
        BXPRC=0.D0;
        BYPRC=0.D0;
        BZPRC=0.D0;
    end
    % C
    if (IOPGEN==0) | (IOPGEN==5),
        HXIMF=0.D0;
        HYIMF=BYIMF;
        HZIMF=BZIMF ;  %! THESE ARE COMPONENTS OF THE PENETRATED FIELD PER UNIT OF THE PENETRATION COEFFICIENT.
        % C                        IN OTHER WORDS, THESE ARE DERIVATIVES OF THE PENETRATION FIELD COMPONENTS WITH RESPECT
        % C                        TO THE PENETRATION COEFFICIENT.   WE ASSUME THAT ONLY THE TRANSVERSE COMPONENT OF THE
        % C                        FIELD PENETRATES INSIDE.
    else
        HXIMF=0.D0;
        HYIMF=0.D0;
        HZIMF=0.D0;
    end
    % C
    % C-----------------------------------------------------------
    % C
    % C    NOW, ADD UP ALL THE COMPONENTS:
    
    DLP1=(PDYN/2.D0)^A(42);
    DLP2=(PDYN/2.D0)^A(43);
    
    TAMP1=A(2)+A(3)*DLP1+A(4)*VBIMF1+A(5)*DST;
    TAMP2=A(6)+A(7)*DLP2+A(8)*VBIMF1+A(9)*DST;
    A_SRC=A(10)+A(11)*DST+A(12)*sqrt(PDYN);
    A_PRC=A(13)+A(14)*DST+A(15)*sqrt(PDYN);
    A_R11=A(16)+A(17)*VBIMF2;
    A_R12=A(18)+A(19)*VBIMF2;
    A_R21=A(20)+A(21)*VBIMF2;
    A_R22=A(22)+A(23)*VBIMF2;
    
    BBX=A(1)*BXCF+TAMP1*BXT1+TAMP2*BXT2+A_SRC*BXSRC+A_PRC*BXPRC ...
        +A_R11*BXR11+A_R12*BXR12+A_R21*BXR21+A_R22*BXR22 ...
        +A(24)*HXIMF+A(25)*HXIMF*STHETAH;
    
    BBY=A(1)*BYCF+TAMP1*BYT1+TAMP2*BYT2+A_SRC*BYSRC+A_PRC*BYPRC ...
        +A_R11*BYR11+A_R12*BYR12+A_R21*BYR21+A_R22*BYR22 ...
        +A(24)*HYIMF+A(25)*HYIMF*STHETAH;
    
    BBZ=A(1)*BZCF+TAMP1*BZT1+TAMP2*BZT2+A_SRC*BZSRC+A_PRC*BZPRC ...
        +A_R11*BZR11+A_R12*BZR12+A_R21*BZR21+A_R22*BZR22 ...
        +A(24)*HZIMF+A(25)*HZIMF*STHETAH;
    
    % C
    % C   AND WE HAVE THE TOTAL EXTERNAL FIELD.
    % C
    % C  NOW, LET US CHECK WHETHER WE HAVE THE CASE (1). IF YES - WE ARE DONE:
    % C
    if (SIGMA<S0-DSIG) ,    %!  (X,Y,Z) IS INSIDE THE MAGNETOSPHERE
        % C-------------------------------------------------------------------------
        BX=BBX;
        BY=BBY;
        BZ=BBZ;
        % C-------------------------------------------------------------------------
    else           %!  THIS IS THE MOST COMPLEX CASE: WE ARE INSIDE
        % C                                             THE INTERPOLATION REGION
        FINT=0.5*(1.-(SIGMA-S0)/DSIG);
        FEXT=0.5*(1.+(SIGMA-S0)/DSIG);
        % C
        [QX,QY,QZ] = T01_DIPOLE (PS,X,Y,Z);
        BX=(BBX+QX)*FINT+OIMFX*FEXT -QX;
        BY=(BBY+QY)*FINT+OIMFY*FEXT -QY;
        BZ=(BBZ+QZ)*FINT+OIMFZ*FEXT -QZ;
        % c
    end  %!   THE CASES (1) AND (2) ARE EXHAUSTED; THE ONLY REMAINING
    % C                      POSSIBILITY IS NOW THE CASE (3):
    % C--------------------------------------------------------------------------
else
    [QX,QY,QZ] = T01_DIPOLE (PS,X,Y,Z);
    BX=OIMFX-QX;
    BY=OIMFY-QY;
    BZ=OIMFZ-QZ;
end % C
% end function EXTALL


% c
% C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
% C
function [BX,BY,BZ] = T01_SHLCAR3X3(X,Y,Z,PS);
%         SUBROUTINE  SHLCAR3X3(X,Y,Z,PS,BX,BY,BZ)
% C
% C   THIS S/R RETURNS THE SHIELDING FIELD FOR THE EARTH'S DIPOLE,
% C   REPRESENTED BY  2x3x3=18 "CARTESIAN" HARMONICS, tilted with respect
% C   to the z=0 plane  (NB#4, p.74)
% C
% C - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% C  The 36 coefficients enter in pairs in the amplitudes of the "cartesian"
% c    harmonics (A(1)-A(36).
% c  The 14 nonlinear parameters (A(37)-A(50) are the scales Pi,Ri,Qi,and Si
% C   entering the arguments of exponents, sines, and cosines in each of the
% C   18 "Cartesian" harmonics  PLUS TWO TILT ANGLES FOR THE CARTESIAN HARMONICS
% C       (ONE FOR THE PSI=0 MODE AND ANOTHER FOR THE PSI=90 MODE)
% C - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% C
%      IMPLICIT  REAL * 8  (A - H, O - Z)
% C
%      DIMENSION A(50)
A = [-901.2327248,895.8011176,817.6208321,-845.5880889, ...
        -83.73539535,86.58542841,336.8781402,-329.3619944,-311.2947120, ...
        308.6011161,31.94469304,-31.30824526,125.8739681,-372.3384278, ...
        -235.4720434,286.7594095,21.86305585,-27.42344605,-150.4874688, ...
        2.669338538,1.395023949,-.5540427503,-56.85224007,3.681827033, ...
        -43.48705106,5.103131905,1.073551279,-.6673083508,12.21404266, ...
        4.177465543,5.799964188,-.3977802319,-1.044652977,.5703560010, ...
        3.536082962,-3.222069852,9.620648151,6.082014949,27.75216226, ...
        12.44199571,5.122226936,6.982039615,20.12149582,6.150973118, ...
        4.663639687,15.73319647,2.303504968,5.840511214,.8385953499E-01, ...
        .3477844929];
% C
P1=A(37);
P2=A(38);
P3=A(39);
R1=A(40);
R2=A(41);
R3=A(42);
Q1=A(43);
Q2=A(44);
Q3=A(45);
S1=A(46);
S2=A(47);
S3=A(48);

T1 =A(49);
T2 =A(50);
% C
CPS=cos(PS);
SPS=sin(PS);
S2PS=2.D0*CPS;      %!   MODIFIED HERE (SIN(2*PS) INSTEAD OF SIN(3*PS))
% C
ST1=sin(PS*T1);
CT1=cos(PS*T1);
ST2=sin(PS*T2);
CT2=cos(PS*T2);

X1=X*CT1-Z*ST1;
Z1=X*ST1+Z*CT1;
X2=X*CT2-Z*ST2;
Z2=X*ST2+Z*CT2;
% C
% C
% c  MAKE THE TERMS IN THE 1ST SUM ("PERPENDICULAR" SYMMETRY):
% C
% C       I=1:
% C
SQPR= sqrt(1.D0/P1^2+1.D0/R1^2);
CYP = cos(Y/P1);
SYP = sin(Y/P1);
CZR = cos(Z1/R1);
SZR = sin(Z1/R1);
EXPR= exp(SQPR*X1);
FX1 =-SQPR*EXPR*CYP*SZR;
HY1 = EXPR/P1*SYP*SZR;
FZ1 =-EXPR*CYP/R1*CZR;
HX1 = FX1*CT1+FZ1*ST1;
HZ1 =-FX1*ST1+FZ1*CT1;

SQPR= sqrt(1.D0/P1^2+1.D0/R2^2);
CYP = cos(Y/P1);
SYP = sin(Y/P1);
CZR = cos(Z1/R2);
SZR = sin(Z1/R2);
EXPR= exp(SQPR*X1);
FX2 =-SQPR*EXPR*CYP*SZR;
HY2 = EXPR/P1*SYP*SZR;
FZ2 =-EXPR*CYP/R2*CZR;
HX2 = FX2*CT1+FZ2*ST1;
HZ2 =-FX2*ST1+FZ2*CT1;

SQPR= sqrt(1.D0/P1^2+1.D0/R3^2);
CYP = cos(Y/P1);
SYP = sin(Y/P1);
CZR = cos(Z1/R3);
SZR = sin(Z1/R3);
EXPR= exp(SQPR*X1);
FX3 =-EXPR*CYP*(SQPR*Z1*CZR+SZR/R3*(X1+1.D0/SQPR));
HY3 = EXPR/P1*SYP*(Z1*CZR+X1/R3*SZR/SQPR);
FZ3 =-EXPR*CYP*(CZR*(1.D0+X1/R3^2/SQPR)-Z1/R3*SZR);
HX3 = FX3*CT1+FZ3*ST1;
HZ3 =-FX3*ST1+FZ3*CT1;
% C
% C       I=2:
% C
SQPR= sqrt(1.D0/P2^2+1.D0/R1^2);
CYP = cos(Y/P2);
SYP = sin(Y/P2);
CZR = cos(Z1/R1);
SZR = sin(Z1/R1);
EXPR= exp(SQPR*X1);
FX4 =-SQPR*EXPR*CYP*SZR;
HY4 = EXPR/P2*SYP*SZR;
FZ4 =-EXPR*CYP/R1*CZR;
HX4 = FX4*CT1+FZ4*ST1;
HZ4 =-FX4*ST1+FZ4*CT1;

SQPR= sqrt(1.D0/P2^2+1.D0/R2^2);
CYP = cos(Y/P2);
SYP = sin(Y/P2);
CZR = cos(Z1/R2);
SZR = sin(Z1/R2);
EXPR= exp(SQPR*X1);
FX5 =-SQPR*EXPR*CYP*SZR;
HY5 = EXPR/P2*SYP*SZR;
FZ5 =-EXPR*CYP/R2*CZR;
HX5 = FX5*CT1+FZ5*ST1;
HZ5 =-FX5*ST1+FZ5*CT1;

SQPR= sqrt(1.D0/P2^2+1.D0/R3^2);
CYP = cos(Y/P2);
SYP = sin(Y/P2);
CZR = cos(Z1/R3);
SZR = sin(Z1/R3);
EXPR= exp(SQPR*X1);
FX6 =-EXPR*CYP*(SQPR*Z1*CZR+SZR/R3*(X1+1.D0/SQPR));
HY6 = EXPR/P2*SYP*(Z1*CZR+X1/R3*SZR/SQPR);
FZ6 =-EXPR*CYP*(CZR*(1.D0+X1/R3^2/SQPR)-Z1/R3*SZR);
HX6 = FX6*CT1+FZ6*ST1;
HZ6 =-FX6*ST1+FZ6*CT1;
% C
% C      I=3:
% C
SQPR= sqrt(1.D0/P3^2+1.D0/R1^2);
CYP = cos(Y/P3);
SYP = sin(Y/P3);
CZR = cos(Z1/R1);
SZR = sin(Z1/R1);
EXPR= exp(SQPR*X1);
FX7 =-SQPR*EXPR*CYP*SZR;
HY7 = EXPR/P3*SYP*SZR;
FZ7 =-EXPR*CYP/R1*CZR;
HX7 = FX7*CT1+FZ7*ST1;
HZ7 =-FX7*ST1+FZ7*CT1;

SQPR= sqrt(1.D0/P3^2+1.D0/R2^2);
CYP = cos(Y/P3);
SYP = sin(Y/P3);
CZR = cos(Z1/R2);
SZR = sin(Z1/R2);
EXPR= exp(SQPR*X1);
FX8 =-SQPR*EXPR*CYP*SZR;
HY8 = EXPR/P3*SYP*SZR;
FZ8 =-EXPR*CYP/R2*CZR;
HX8 = FX8*CT1+FZ8*ST1;
HZ8 =-FX8*ST1+FZ8*CT1;

SQPR= sqrt(1.D0/P3^2+1.D0/R3^2);
CYP = cos(Y/P3);
SYP = sin(Y/P3);
CZR = cos(Z1/R3);
SZR = sin(Z1/R3);
EXPR= exp(SQPR*X1);
FX9 =-EXPR*CYP*(SQPR*Z1*CZR+SZR/R3*(X1+1.D0/SQPR));
HY9 = EXPR/P3*SYP*(Z1*CZR+X1/R3*SZR/SQPR);
FZ9 =-EXPR*CYP*(CZR*(1.D0+X1/R3^2/SQPR)-Z1/R3*SZR);
HX9 = FX9*CT1+FZ9*ST1;
HZ9 =-FX9*ST1+FZ9*CT1;
% C
A1=A(1)+A(2)*CPS;
A2=A(3)+A(4)*CPS;
A3=A(5)+A(6)*CPS;
A4=A(7)+A(8)*CPS;
A5=A(9)+A(10)*CPS;
A6=A(11)+A(12)*CPS;
A7=A(13)+A(14)*CPS;
A8=A(15)+A(16)*CPS;
A9=A(17)+A(18)*CPS;
BX=A1*HX1+A2*HX2+A3*HX3+A4*HX4+A5*HX5+A6*HX6+A7*HX7+A8*HX8+A9*HX9;
BY=A1*HY1+A2*HY2+A3*HY3+A4*HY4+A5*HY5+A6*HY6+A7*HY7+A8*HY8+A9*HY9;
BZ=A1*HZ1+A2*HZ2+A3*HZ3+A4*HZ4+A5*HZ5+A6*HZ6+A7*HZ7+A8*HZ8+A9*HZ9;


% c  MAKE THE TERMS IN THE 2ND SUM ("PARALLEL" SYMMETRY):
% C
% C       I=1
% C
SQQS= sqrt(1.D0/Q1^2+1.D0/S1^2);
CYQ = cos(Y/Q1);
SYQ = sin(Y/Q1);
CZS = cos(Z2/S1);
SZS = sin(Z2/S1);
EXQS= exp(SQQS*X2);
FX1 =-SQQS*EXQS*CYQ*CZS *SPS;
HY1 = EXQS/Q1*SYQ*CZS   *SPS;
FZ1 = EXQS*CYQ/S1*SZS   *SPS;
HX1 = FX1*CT2+FZ1*ST2;
HZ1 =-FX1*ST2+FZ1*CT2;

SQQS= sqrt(1.D0/Q1^2+1.D0/S2^2);
CYQ = cos(Y/Q1);
SYQ = sin(Y/Q1);
CZS = cos(Z2/S2);
SZS = sin(Z2/S2);
EXQS= exp(SQQS*X2);
FX2 =-SQQS*EXQS*CYQ*CZS *SPS;
HY2 = EXQS/Q1*SYQ*CZS   *SPS;
FZ2 = EXQS*CYQ/S2*SZS   *SPS;
HX2 = FX2*CT2+FZ2*ST2;
HZ2 =-FX2*ST2+FZ2*CT2;

SQQS= sqrt(1.D0/Q1^2+1.D0/S3^2);
CYQ = cos(Y/Q1);
SYQ = sin(Y/Q1);
CZS = cos(Z2/S3);
SZS = sin(Z2/S3);
EXQS= exp(SQQS*X2);
FX3 =-SQQS*EXQS*CYQ*CZS *SPS;
HY3 = EXQS/Q1*SYQ*CZS   *SPS;
FZ3 = EXQS*CYQ/S3*SZS   *SPS;
HX3 = FX3*CT2+FZ3*ST2;
HZ3 =-FX3*ST2+FZ3*CT2;
% C
% C       I=2:
% C
SQQS= sqrt(1.D0/Q2^2+1.D0/S1^2);
CYQ = cos(Y/Q2);
SYQ = sin(Y/Q2);
CZS = cos(Z2/S1);
SZS = sin(Z2/S1);
EXQS= exp(SQQS*X2);
FX4 =-SQQS*EXQS*CYQ*CZS *SPS;
HY4 = EXQS/Q2*SYQ*CZS   *SPS;
FZ4 = EXQS*CYQ/S1*SZS   *SPS;
HX4 = FX4*CT2+FZ4*ST2;
HZ4 =-FX4*ST2+FZ4*CT2;

SQQS= sqrt(1.D0/Q2^2+1.D0/S2^2);
CYQ = cos(Y/Q2);
SYQ = sin(Y/Q2);
CZS = cos(Z2/S2);
SZS = sin(Z2/S2);
EXQS= exp(SQQS*X2);
FX5 =-SQQS*EXQS*CYQ*CZS *SPS;
HY5 = EXQS/Q2*SYQ*CZS   *SPS;
FZ5 = EXQS*CYQ/S2*SZS   *SPS;
HX5 = FX5*CT2+FZ5*ST2;
HZ5 =-FX5*ST2+FZ5*CT2;

SQQS= sqrt(1.D0/Q2^2+1.D0/S3^2);
CYQ = cos(Y/Q2);
SYQ = sin(Y/Q2);
CZS = cos(Z2/S3);
SZS = sin(Z2/S3);
EXQS= exp(SQQS*X2);
FX6 =-SQQS*EXQS*CYQ*CZS *SPS;
HY6 = EXQS/Q2*SYQ*CZS   *SPS;
FZ6 = EXQS*CYQ/S3*SZS   *SPS;
HX6 = FX6*CT2+FZ6*ST2;
HZ6 =-FX6*ST2+FZ6*CT2;
% C
% C       I=3:
% C
SQQS= sqrt(1.D0/Q3^2+1.D0/S1^2);
CYQ = cos(Y/Q3);
SYQ = sin(Y/Q3);
CZS = cos(Z2/S1);
SZS = sin(Z2/S1);
EXQS= exp(SQQS*X2);
FX7 =-SQQS*EXQS*CYQ*CZS *SPS;
HY7 = EXQS/Q3*SYQ*CZS   *SPS;
FZ7 = EXQS*CYQ/S1*SZS   *SPS;
HX7 = FX7*CT2+FZ7*ST2;
HZ7 =-FX7*ST2+FZ7*CT2;

SQQS= sqrt(1.D0/Q3^2+1.D0/S2^2);
CYQ = cos(Y/Q3);
SYQ = sin(Y/Q3);
CZS = cos(Z2/S2);
SZS = sin(Z2/S2);
EXQS= exp(SQQS*X2);
FX8 =-SQQS*EXQS*CYQ*CZS *SPS;
HY8 = EXQS/Q3*SYQ*CZS   *SPS;
FZ8 = EXQS*CYQ/S2*SZS   *SPS;
HX8 = FX8*CT2+FZ8*ST2;
HZ8 =-FX8*ST2+FZ8*CT2;

SQQS= sqrt(1.D0/Q3^2+1.D0/S3^2);
CYQ = cos(Y/Q3);
SYQ = sin(Y/Q3);
CZS = cos(Z2/S3);
SZS = sin(Z2/S3);
EXQS= exp(SQQS*X2);
FX9 =-SQQS*EXQS*CYQ*CZS *SPS;
HY9 = EXQS/Q3*SYQ*CZS   *SPS;
FZ9 = EXQS*CYQ/S3*SZS   *SPS;
HX9 = FX9*CT2+FZ9*ST2;
HZ9 =-FX9*ST2+FZ9*CT2;

A1=A(19)+A(20)*S2PS;
A2=A(21)+A(22)*S2PS;
A3=A(23)+A(24)*S2PS;
A4=A(25)+A(26)*S2PS;
A5=A(27)+A(28)*S2PS;
A6=A(29)+A(30)*S2PS;
A7=A(31)+A(32)*S2PS;
A8=A(33)+A(34)*S2PS;
A9=A(35)+A(36)*S2PS;

BX=BX+A1*HX1+A2*HX2+A3*HX3+A4*HX4+A5*HX5+A6*HX6+A7*HX7+A8*HX8+A9*HX9;
BY=BY+A1*HY1+A2*HY2+A3*HY3+A4*HY4+A5*HY5+A6*HY6+A7*HY7+A8*HY8+A9*HY9;
BZ=BZ+A1*HZ1+A2*HZ2+A3*HZ3+A4*HZ4+A5*HZ5+A6*HZ6+A7*HZ7+A8*HZ8+A9*HZ9;
% C
% end of function SHLCAR3x3


% c
% c############################################################################
% c
% C
function [BX1,BY1,BZ1,BX2,BY2,BZ2] = T01_DEFORMED(IOPT,PS,X,Y,Z)
%      SUBROUTINE DEFORMED (IOPT,PS,X,Y,Z,BX1,BY1,BZ1,BX2,BY2,BZ2)
% C
% C   IOPT - TAIL FIELD MODE FLAG:   IOPT=0 - THE TWO TAIL MODES ARE ADDED UP
% C                                  IOPT=1 - MODE 1 ONLY
% C                                  IOPT=2 - MODE 2 ONLY
% C
% C   CALCULATES GSM COMPONENTS OF TWO UNIT-AMPLITUDE TAIL FIELD MODES,
% C    TAKING INTO ACCOUNT BOTH EFFECTS OF DIPOLE TILT:
% C    WARPING IN Y-Z (DONE BY THE S/R WARPED) AND BENDING IN X-Z (DONE BY THIS SUBROUTINE)
% C
%      IMPLICIT REAL*8 (A-H,O-Z)   
%      COMMON /RH0/ RH0
global T01;
%      DATA RH2,IEPS /-5.2D0,3/
RH2 = -5.2D0;
IEPS=3;
% C
% C  RH0,RH1,RH2, AND IEPS CONTROL THE TILT-RELATED DEFORMATION OF THE TAIL FIELD
% C
SPS=sin(PS);
CPS=sqrt(1.D0-SPS^2);
R2=X^2+Y^2+Z^2;
R=sqrt(R2);
ZR=Z/R;
RH=T01.RH0+RH2*ZR^2;
DRHDR=-ZR/R*2.D0*RH2*ZR;
DRHDZ= 2.D0*RH2*ZR/R;
% C
RRH=R/RH;
F=1.D0/(1.D0+RRH^IEPS)^(1.D0/IEPS);
DFDR=-RRH^(IEPS-1)*F^(IEPS+1)/RH;
DFDRH=-RRH*DFDR;
% c
SPSAS=SPS*F;
CPSAS=sqrt(1.D0-SPSAS^2);
% C
XAS=X*CPSAS-Z*SPSAS;
ZAS=X*SPSAS+Z*CPSAS;
% C
FACPS=SPS/CPSAS*(DFDR+DFDRH*DRHDR)/R;
PSASX=FACPS*X;
PSASY=FACPS*Y;
PSASZ=FACPS*Z+SPS/CPSAS*DFDRH*DRHDZ;
% C
DXASDX=CPSAS-ZAS*PSASX;
DXASDY=-ZAS*PSASY;
DXASDZ=-SPSAS-ZAS*PSASZ;
DZASDX=SPSAS+XAS*PSASX;
DZASDY=XAS*PSASY;
DZASDZ=CPSAS+XAS*PSASZ;
FAC1=DXASDZ*DZASDY-DXASDY*DZASDZ;
FAC2=DXASDX*DZASDZ-DXASDZ*DZASDX;
FAC3=DZASDX*DXASDY-DXASDX*DZASDY;
% C
% C     DEFORM:
% C
[BXAS1,BYAS1,BZAS1,BXAS2,BYAS2,BZAS2] = T01_WARPED(IOPT,PS,XAS,Y,ZAS);
% C
BX1=BXAS1*DZASDZ-BZAS1*DXASDZ +BYAS1*FAC1;
BY1=BYAS1*FAC2;
BZ1=BZAS1*DXASDX-BXAS1*DZASDX +BYAS1*FAC3;

BX2=BXAS2*DZASDZ-BZAS2*DXASDZ +BYAS2*FAC1;
BY2=BYAS2*FAC2;
BZ2=BZAS2*DXASDX-BXAS2*DZASDX +BYAS2*FAC3;

% end of function DEFORMED


% C
% C------------------------------------------------------------------
% C
function [BX1,BY1,BZ1,BX2,BY2,BZ2] = T01_WARPED(IOPT,PS,X,Y,Z)
%      SUBROUTINE WARPED (IOPT,PS,X,Y,Z,BX1,BY1,BZ1,BX2,BY2,BZ2)
% C
% C   CALCULATES GSM COMPONENTS OF THE WARPED FIELD FOR TWO TAIL UNIT MODES.
% C   THE WARPING DEFORMATION IS IMPOSED ON THE UNWARPED FIELD, COMPUTED
% C   BY THE S/R "UNWARPED".  THE WARPING PARAMETER G WAS OBTAINED BY LEAST
% C   SQUARES FITTING TO THE ENTIRE DATASET.
% C
% C   IOPT - TAIL FIELD MODE FLAG:   IOPT=0 - THE TWO TAIL MODES ARE ADDED UP
% C                                  IOPT=1 - MODE 1 ONLY
% C                                  IOPT=2 - MODE 2 ONLY
% C
%      IMPLICIT REAL*8 (A-H,O-Z)
% C
%      COMMON /G/ G
global T01;
DGDX=0.D0;
XL=20.D0;
DXLDX=0.D0;

SPS=sin(PS);
RHO2=Y^2+Z^2;
RHO=sqrt(RHO2);

if (Y==0.D0) & (Z==0.D0),
    PHI=0.D0;
    CPHI=1.D0;
    SPHI=0.D0;
else
    PHI=atan2(Z,Y);
    CPHI=Y/RHO;
    SPHI=Z/RHO;
end

RR4L4=RHO/(RHO2^2+XL^4);

F=PHI+T01.G*RHO2*RR4L4*CPHI*SPS;
DFDPHI=1.D0-T01.G*RHO2*RR4L4*SPHI*SPS;
DFDRHO=T01.G*RR4L4^2*(3.D0*XL^4-RHO2^2)*CPHI*SPS;
DFDX=RR4L4*CPHI*SPS*(DGDX*RHO2-T01.G*RHO*RR4L4*4.D0*XL^3*DXLDX);

CF=cos(F);
SF=sin(F);
YAS=RHO*CF;
ZAS=RHO*SF;

[BX_AS1,BY_AS1,BZ_AS1,BX_AS2,BY_AS2,BZ_AS2] = T01_UNWARPED(IOPT,X,YAS,ZAS);

BRHO_AS =  BY_AS1*CF+BZ_AS1*SF;      %!   DEFORM THE 1ST MODE
BPHI_AS = -BY_AS1*SF+BZ_AS1*CF;

BRHO_S = BRHO_AS*DFDPHI;
BPHI_S = BPHI_AS-RHO*(BX_AS1*DFDX+BRHO_AS*DFDRHO);
BX1    = BX_AS1*DFDPHI;

BY1    = BRHO_S*CPHI-BPHI_S*SPHI;
BZ1    = BRHO_S*SPHI+BPHI_S*CPHI;    %!   DONE

BRHO_AS =  BY_AS2*CF+BZ_AS2*SF;      %!   DEFORM THE 2ND MODE
BPHI_AS = -BY_AS2*SF+BZ_AS2*CF;

BRHO_S = BRHO_AS*DFDPHI;
BPHI_S = BPHI_AS-RHO*(BX_AS2*DFDX+BRHO_AS*DFDRHO);
BX2    = BX_AS2*DFDPHI;

BY2    = BRHO_S*CPHI-BPHI_S*SPHI;
BZ2    = BRHO_S*SPHI+BPHI_S*CPHI;    %!   DONE

% end of function WARPED

% C
% C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% C
function [BX1,BY1,BZ1,BX2,BY2,BZ2] = T01_UNWARPED(IOPT,X,Y,Z);
%       SUBROUTINE UNWARPED (IOPT,X,Y,Z,BX1,BY1,BZ1,BX2,BY2,BZ2)

% C   IOPT - TAIL FIELD MODE FLAG:   IOPT=0 - THE TWO TAIL MODES ARE ADDED UP
% C                                  IOPT=1 - MODE 1 ONLY
% C                                  IOPT=2 - MODE 2 ONLY
% C
% C    CALCULATES GSM COMPONENTS OF THE SHIELDED FIELD OF TWO TAIL MODES WITH UNIT
% C    AMPLITUDES,  WITHOUT ANY WARPING OR BENDING.  NONLINEAR PARAMETERS OF THE MODES
% C    ARE FORWARDED HERE VIA A COMMON BLOCK /TAIL/.
% C
%      IMPLICIT REAL*8 (A-H,O-Z)
% C
%      DIMENSION A1(60),A2(60)  %!   TAIL SHIELDING FIELD PARAMETERS FOR THE MODES #1 & #2

%      COMMON /TAIL/ DXSHIFT1,DXSHIFT2,D0,DELTADY
global T01;
% D0 -> T0.TAIL.D
% C
%      DATA DELTADX1,ALPHA1,XSHIFT1 /1.D0,1.1D0,6.D0/
DELTADX1 = 1.D0;
ALPHA1 = 1.1D0;
XSHIFT1 = 6.D0;
%      DATA DELTADX2,ALPHA2,XSHIFT2 /0.D0,.25D0,4.D0/
DELTADX2 = 0.D0;
ALPHA2 = .25D0;
XSHIFT2 = 4.D0;

A1 = [-25.45869857,57.35899080,317.5501869,-2.626756717, ...
        -93.38053698,-199.6467926,-858.8129729,34.09192395,845.4214929, ...
        -29.07463068,47.10678547,-128.9797943,-781.7512093,6.165038619, ...
        167.8905046,492.0680410,1654.724031,-46.77337920,-1635.922669, ...
        40.86186772,-.1349775602,-.9661991179E-01,-.1662302354, ...
        .002810467517,.2487355077,.1025565237,-14.41750229,-.8185333989, ...
        11.07693629,.7569503173,-9.655264745,112.2446542,777.5948964, ...
        -5.745008536,-83.03921993,-490.2278695,-1155.004209,39.08023320, ...
        1172.780574,-39.44349797,-14.07211198,-40.41201127,-313.2277343, ...
        2.203920979,8.232835341,197.7065115,391.2733948,-18.57424451, ...
        -437.2779053,23.04976898,11.75673963,13.60497313,4.691927060, ...
        18.20923547,27.59044809,6.677425469,1.398283308,2.839005878, ...
        31.24817706,24.53577264];

A2 = [-287187.1962,4970.499233,410490.1952,-1347.839052, ...
        -386370.3240,3317.983750,-143462.3895,5706.513767,171176.2904, ...
        250.8882750,-506570.8891,5733.592632,397975.5842,9771.762168, ...
        -941834.2436,7990.975260,54313.10318,447.5388060,528046.3449, ...
        12751.04453,-21920.98301,-21.05075617,31971.07875,3012.641612, ...
        -301822.9103,-3601.107387,1797.577552,-6.315855803,142578.8406, ...
        13161.93640,804184.8410,-14168.99698,-851926.6360,-1890.885671, ...
        972475.6869,-8571.862853,26432.49197,-2554.752298,-482308.3431, ...
        -4391.473324,105155.9160,-1134.622050,-74353.53091,-5382.670711, ...
        695055.0788,-916.3365144,-12111.06667,67.20923358,-367200.9285, ...
        -21414.14421,14.75567902,20.75638190,59.78601609,16.86431444, ...
        32.58482365,23.69472951,17.24977936,13.64902647,68.40989058, ...
        11.67828167];

%      DATA XM1,XM2/2*-12.D0/
XM1 = -12;
XM2 = -12;

if (IOPT~=2),
    
    XSC1=(X-XSHIFT1-T01.TAIL.DXSHIFT1)*ALPHA1-XM1*(ALPHA1-1.D0);
    YSC1=Y*ALPHA1;
    ZSC1=Z*ALPHA1;
    D0SC1=T01.TAIL.D*ALPHA1;   %! HERE WE USE A SINGLE VALUE D0 OF THE THICKNESS FOR BOTH MODES
    
    [FX1,FY1,FZ1] = T01_TAILDISK(D0SC1,DELTADX1,T01.TAIL.DELTADY,XSC1,YSC1,ZSC1);
    [HX1,HY1,HZ1] = T01_SHLCAR5X5(A1,X,Y,Z,T01.TAIL.DXSHIFT1);
    
    BX1=FX1+HX1;
    BY1=FY1+HY1;
    BZ1=FZ1+HZ1;
    
    if (IOPT==1),
        BX2=0.D0;
        BY2=0.D0;
        BZ2=0.D0;
        return
    end
    
end

XSC2=(X-XSHIFT2-T01.TAIL.DXSHIFT2)*ALPHA2-XM2*(ALPHA2-1.D0);
YSC2=Y*ALPHA2;
ZSC2=Z*ALPHA2;
D0SC2=T01.TAIL.D*ALPHA2;   %! HERE WE USE A SINGLE VALUE D0 OF THE THICKNESS FOR BOTH MODES

[FX2,FY2,FZ2] = T01_TAILDISK(D0SC2,DELTADX2,T01.TAIL.DELTADY,XSC2,YSC2,ZSC2);
[HX2,HY2,HZ2] = T01_SHLCAR5X5(A2,X,Y,Z,T01.TAIL.DXSHIFT2);

BX2=FX2+HX2;
BY2=FY2+HY2;
BZ2=FZ2+HZ2;

if (IOPT==2),
    BX1=0.D0;
    BY1=0.D0;
    BZ1=0.D0;
    return
end

% end of function UNWARPED

% C
% C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
% C
function [BX,BY,BZ] = T01_TAILDISK(D0,DELTADX,DELTADY,X,Y,Z);
%      SUBROUTINE TAILDISK(D0,DELTADX,DELTADY,X,Y,Z,BX,BY,BZ)
% c
% c       THIS SUBROUTINE COMPUTES THE COMPONENTS OF THE TAIL CURRENT FIELD,
% C       SIMILAR TO THAT DESCRIBED BY TSYGANENKO AND PEREDO (1994).  THE
% C       DIFFERENCE IS THAT NOW WE USE SPACEWARPING, AS DESCRIBED IN OUR
% C       PAPER ON MODELING BIRKELAND CURRENTS (TSYGANENKO AND STERN, 1996)
% C       INSTEAD OF SHEARING IT IN THE SPIRIT OF THE T89 TAIL MODEL.
% C
%      IMPLICIT REAL*8 (A-H,O-Z)
% c
%      DIMENSION F(5),B(5),C(5)
% C
F = [-71.09346626D0,-1014.308601D0,-1272.939359D0, ...
        -3224.935936D0,-44546.86232D0];
B = [10.90101242D0,12.68393898D0,13.51791954D0,14.86775017D0, ...
        15.12306404D0];
C = [ .7954069972D0,.6716601849D0,1.174866319D0,2.565249920D0, ...
        10.01986790D0];
% C
RHO=sqrt(X^2+Y^2);
DRHODX=X/RHO;
DRHODY=Y/RHO;

DEX=exp(X/7.D0);
D=D0+DELTADY*(Y/20.D0)^2  +DELTADX*DEX; %!   THE LAST TERM (INTRODUCED 10/11/2000) MAKES THE SHEET
DDDY=DELTADY*Y*0.005D0;                  %!   THICKEN SUNWARD, TO AVOID PROBLEMS IN THE SUBSOLAR REGION
DDDX=DELTADX/7.D0*DEX;
% C
DZETA=sqrt(Z^2+D^2);  %!  THIS IS THE SAME SIMPLE WAY TO SPREAD
% C                                        OUT THE SHEET, AS THAT USED IN T89
DDZETADX=D*DDDX/DZETA;
DDZETADY=D*DDDY/DZETA;
DDZETADZ=Z/DZETA;

% C
DBX=0.D0;
DBY=0.D0;
DBZ=0.D0;
% C
for I=1:5,
    % C
    BI=B(I);
    CI=C(I);
    % C
    S1=sqrt((RHO+BI)^2+(DZETA+CI)^2);
    S2=sqrt((RHO-BI)^2+(DZETA+CI)^2);
    
    DS1DRHO=(RHO+BI)/S1;
    DS2DRHO=(RHO-BI)/S2;
    DS1DDZ=(DZETA+CI)/S1;
    DS2DDZ=(DZETA+CI)/S2;
    % C
    DS1DX=DS1DRHO*DRHODX  +DS1DDZ*DDZETADX;
    DS1DY=DS1DRHO*DRHODY  +   DS1DDZ*DDZETADY;
    DS1DZ=                      DS1DDZ*DDZETADZ;
    % C
    DS2DX=DS2DRHO*DRHODX  +DS2DDZ*DDZETADX;
    DS2DY=DS2DRHO*DRHODY  +   DS2DDZ*DDZETADY;
    DS2DZ=                    DS2DDZ*DDZETADZ;
    % C
    S1TS2=S1*S2;
    S1PS2=S1+S2;
    S1PS2SQ=S1PS2^2;
    
    FAC1=sqrt(S1PS2SQ-(2.D0*BI)^2);
    AS=FAC1/(S1TS2*S1PS2SQ);
    DASDS1=(1.D0/(FAC1*S2)-AS/S1PS2*(S2*S2+S1*(3.D0*S1+4.D0*S2))) ...
        /(S1*S1PS2);
    DASDS2=(1.D0/(FAC1*S1)-AS/S1PS2*(S1*S1+S2*(3.D0*S2+4.D0*S1))) ...
        /(S2*S1PS2);
    % C
    DASDX=DASDS1*DS1DX+DASDS2*DS2DX;
    DASDY=DASDS1*DS1DY+DASDS2*DS2DY;
    DASDZ=DASDS1*DS1DZ+DASDS2*DS2DZ;
    % C
    DBX=DBX-F(I)*X*DASDZ;
    DBY=DBY-F(I)*Y*DASDZ;
    DBZ=DBZ+F(I)*(2.D0*AS+X*DASDX+Y*DASDY);
end

BX=DBX;
BY=DBY;
BZ=DBZ;

% end of function TAILDISK


% C
% C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
% C
% C THIS CODE RETURNS THE SHIELDING FIELD REPRESENTED BY  5x5=25 "CARTESIAN"
% C    HARMONICS
% C
function [HX,HY,HZ] = T01_SHLCAR5X5(A,X,Y,Z,DSHIFT);
%         SUBROUTINE  SHLCAR5X5(A,X,Y,Z,DSHIFT,HX,HY,HZ)
% C
% C - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% C  The NLIN coefficients are the amplitudes of the "cartesian"
% c    harmonics (A(1)-A(NLIN).
% c  The NNP nonlinear parameters (A(NLIN+1)-A(NTOT) are the scales Pi and Ri
% C   entering the arguments of exponents, sines, and cosines in each of the
% C   NLIN "Cartesian" harmonics
% C - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% C
%         IMPLICIT  REAL * 8  (A - H, O - Z)
% C
%         DIMENSION A(60)
% C
DHX=0.D0;
DHY=0.D0;
DHZ=0.D0;

L=0;
% C
for I=1:5,
    %         DO 2 I=1,5
    RP=1.D0/A(50+I);
    CYPI=cos(Y*RP);
    SYPI=sin(Y*RP);
    % C
    for K=1:5,
        %         DO 2 K=1,5
        RR=1.D0/A(55+K);
        SZRK=sin(Z*RR);
        CZRK=cos(Z*RR);
        SQPR=sqrt(RP^2+RR^2);
        EPR=exp(X*SQPR);
        % C
        DBX=-SQPR*EPR*CYPI*SZRK;
        DBY= RP*EPR*SYPI*SZRK;
        DBZ=-RR*EPR*CYPI*CZRK;
        
        L=L+2;
        COEF=A(L-1)+A(L)*DSHIFT;
        
        DHX=DHX+COEF*DBX;
        DHY=DHY+COEF*DBY;
        DHZ=DHZ+COEF*DBZ;
        % c
    end %  2      CONTINUE
end 

HX=DHX;
HY=DHY;
HZ=DHZ;
% C

% end of function SHLCAR5X5


% c
% c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% C
function [BX11,BY11,BZ11,BX12,BY12,BZ12,...
        BX21,BY21,BZ21,BX22,BY22,BZ22] = T01_BIRK_TOT (IOPB,PS,X,Y,Z);
%      SUBROUTINE BIRK_TOT (IOPB,PS,X,Y,Z,BX11,BY11,BZ11,BX12,BY12,BZ12,
%     *                          BX21,BY21,BZ21,BX22,BY22,BZ22)
% C
% C      IOPB -  BIRKELAND FIELD MODE FLAG:
% C         IOPB=0 - ALL COMPONENTS
% C         IOPB=1 - REGION 1, MODES 1 & 2
% C         IOPB=2 - REGION 2, MODES 1 & 2
% C
%      IMPLICIT REAL*8 (A-H,O-Z)
%      DIMENSION SH11(86),SH12(86),SH21(86),SH22(86)
global T01;    
%      COMMON /BIRKPAR/ XKAPPA1,XKAPPA2   %!  INPUT PARAMETERS, SPECIFIED FROM S/R EXTALL
%      COMMON /DPHI_B_RHO0/ DPHI,B,RHO_0,XKAPPA %! PARAMETERS, CONTROLLING THE DAY-NIGHT ASYMMETRY OF F.A.C.

SH11 = [46488.84663,-15541.95244,-23210.09824,-32625.03856, ...
        -109894.4551,-71415.32808,58168.94612,55564.87578,-22890.60626, ...
        -6056.763968,5091.368100,239.7001538,-13899.49253,4648.016991, ...
        6971.310672,9699.351891,32633.34599,21028.48811,-17395.96190, ...
        -16461.11037,7447.621471,2528.844345,-1934.094784,-588.3108359, ...
        -32588.88216,10894.11453,16238.25044,22925.60557,77251.11274, ...
        50375.97787,-40763.78048,-39088.60660,15546.53559,3559.617561, ...
        -3187.730438,309.1487975,88.22153914,-243.0721938,-63.63543051, ...
        191.1109142,69.94451996,-187.9539415,-49.89923833,104.0902848, ...
        -120.2459738,253.5572433,89.25456949,-205.6516252,-44.93654156, ...
        124.7026309,32.53005523,-98.85321751,-36.51904756,98.88241690, ...
        24.88493459,-55.04058524,61.14493565,-128.4224895,-45.35023460, ...
        105.0548704,-43.66748755,119.3284161,31.38442798,-92.87946767, ...
        -33.52716686,89.98992001,25.87341323,-48.86305045,59.69362881, ...
        -126.5353789,-44.39474251,101.5196856,59.41537992,41.18892281, ...
        80.86101200,3.066809418,7.893523804,30.56212082,10.36861082, ...
        8.222335945,19.97575641,2.050148531,4.992657093,2.300564232, ...
        .2256245602,-.05841594319];

SH12 = [210260.4816,-1443587.401,-1468919.281,281939.2993, ...
        -1131124.839,729331.7943,2573541.307,304616.7457,468887.5847, ...
        181554.7517,-1300722.650,-257012.8601,645888.8041,-2048126.412, ...
        -2529093.041,571093.7972,-2115508.353,1122035.951,4489168.802, ...
        75234.22743,823905.6909,147926.6121,-2276322.876,-155528.5992, ...
        -858076.2979,3474422.388,3986279.931,-834613.9747,3250625.781, ...
        -1818680.377,-7040468.986,-414359.6073,-1295117.666,-346320.6487, ...
        3565527.409,430091.9496,-.1565573462,7.377619826,.4115646037, ...
        -6.146078880,3.808028815,-.5232034932,1.454841807,-12.32274869, ...
        -4.466974237,-2.941184626,-.6172620658,12.64613490,1.494922012, ...
        -21.35489898,-1.652256960,16.81799898,-1.404079922,-24.09369677, ...
        -10.99900839,45.94237820,2.248579894,31.91234041,7.575026816, ...
        -45.80833339,-1.507664976,14.60016998,1.348516288,-11.05980247, ...
        -5.402866968,31.69094514,12.28261196,-37.55354174,4.155626879, ...
        -33.70159657,-8.437907434,36.22672602,145.0262164,70.73187036, ...
        85.51110098,21.47490989,24.34554406,31.34405345,4.655207476, ...
        5.747889264,7.802304187,1.844169801,4.867254550,2.941393119, ...
        .1379899178,.06607020029];

SH21 = [162294.6224,503885.1125,-27057.67122,-531450.1339, ...
        84747.05678,-237142.1712,84133.61490,259530.0402,69196.05160, ...
        -189093.5264,-19278.55134,195724.5034,-263082.6367,-818899.6923, ...
        43061.10073,863506.6932,-139707.9428,389984.8850,-135167.5555, ...
        -426286.9206,-109504.0387,295258.3531,30415.07087,-305502.9405, ...
        100785.3400,315010.9567,-15999.50673,-332052.2548,54964.34639, ...
        -152808.3750,51024.67566,166720.0603,40389.67945,-106257.7272, ...
        -11126.14442,109876.2047,2.978695024,558.6019011,2.685592939, ...
        -338.0004730,-81.99724090,-444.1102659,89.44617716,212.0849592, ...
        -32.58562625,-982.7336105,-35.10860935,567.8931751,-1.917212423, ...
        -260.2023543,-1.023821735,157.5533477,23.00200055,232.0603673, ...
        -36.79100036,-111.9110936,18.05429984,447.0481000,15.10187415, ...
        -258.7297813,-1.032340149,-298.6402478,-1.676201415,180.5856487, ...
        64.52313024,209.0160857,-53.85574010,-98.52164290,14.35891214, ...
        536.7666279,20.09318806,-309.7349530,58.54144539,67.45226850, ...
        97.92374406,4.752449760,10.46824379,32.91856110,12.05124381, ...
        9.962933904,15.91258637,1.804233877,6.578149088,2.515223491, ...
        .1930034238,-.02261109942];

SH22 = [-131287.8986,-631927.6885,-318797.4173,616785.8782, ...
        -50027.36189,863099.9833,47680.20240,-1053367.944,-501120.3811, ...
        -174400.9476,222328.6873,333551.7374,-389338.7841,-1995527.467, ...
        -982971.3024,1960434.268,297239.7137,2676525.168,-147113.4775, ...
        -3358059.979,-2106979.191,-462827.1322,1017607.960,1039018.475, ...
        520266.9296,2627427.473,1301981.763,-2577171.706,-238071.9956, ...
        -3539781.111,94628.16420,4411304.724,2598205.733,637504.9351, ...
        -1234794.298,-1372562.403,-2.646186796,-31.10055575,2.295799273, ...
        19.20203279,30.01931202,-302.1028550,-14.78310655,162.1561899, ...
        .4943938056,176.8089129,-.2444921680,-100.6148929,9.172262228, ...
        137.4303440,-8.451613443,-84.20684224,-167.3354083,1321.830393, ...
        76.89928813,-705.7586223,18.28186732,-770.1665162,-9.084224422, ...
        436.3368157,-6.374255638,-107.2730177,6.080451222,65.53843753, ...
        143.2872994,-1028.009017,-64.22739330,547.8536586,-20.58928632, ...
        597.3893669,10.17964133,-337.7800252,159.3532209,76.34445954, ...
        84.74398828,12.76722651,27.63870691,32.69873634,5.145153451, ...
        6.310949163,6.996159733,1.971629939,4.436299219,2.904964304, ...
        .1486276863,.06859991529];
% C
T01.DPHI_B_RHO0.XKAPPA=T01.BIRKPAR.XKAPPA1;        %!  FORWARDED IN BIRK_1N2
X_SC=T01.BIRKPAR.XKAPPA1-1.1D0;    %!  FORWARDED IN BIRK_SHL

if (IOPB==0) | (IOPB==1),
    
    [FX11,FY11,FZ11] = T01_BIRK_1N2(1,1,PS,X,Y,Z);           %!  REGION 1, MODE 1
    [HX11,HY11,HZ11] = T01_BIRK_SHL(SH11,PS,X_SC,X,Y,Z);
    BX11=FX11+HX11;
    BY11=FY11+HY11;
    BZ11=FZ11+HZ11;
    
    [FX12,FY12,FZ12] = T01_BIRK_1N2 (1,2,PS,X,Y,Z);           %!  REGION 1, MODE 2
    [HX12,HY12,HZ12] = T01_BIRK_SHL (SH12,PS,X_SC,X,Y,Z);
    BX12=FX12+HX12;
    BY12=FY12+HY12;
    BZ12=FZ12+HZ12;
    
end

T01.DPHI_B_RHO0.XKAPPA=T01.BIRKPAR.XKAPPA2  ;      %!  FORWARDED IN BIRK_1N2
X_SC=T01.BIRKPAR.XKAPPA2-1.0D0 ;   %!  FORWARDED IN BIRK_SHL

if (IOPB==0) | (IOPB==2),
    
    [FX21,FY21,FZ21] = T01_BIRK_1N2(2,1,PS,X,Y,Z);           %!  REGION 2, MODE 1
    [HX21,HY21,HZ21] = T01_BIRK_SHL(SH21,PS,X_SC,X,Y,Z);
    BX21=FX21+HX21;
    BY21=FY21+HY21;
    BZ21=FZ21+HZ21;
    
    [FX22,FY22,FZ22] = T01_BIRK_1N2(2,2,PS,X,Y,Z);           %!  REGION 2, MODE 2
    [HX22,HY22,HZ22] = T01_BIRK_SHL(SH22,PS,X_SC,X,Y,Z);
    BX22=FX22+HX22;
    BY22=FY22+HY22;
    BZ22=FZ22+HZ22;
    
end

% end of function BIRK_TOT

% C
% c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% c
% c
function [BX,BY,BZ] = T01_BIRK_1N2 (NUMB,MODE,PS,X,Y,Z);        %!   NB# 6, P.60
%      SUBROUTINE BIRK_1N2 (NUMB,MODE,PS,X,Y,Z,BX,BY,BZ)        %!   NB# 6, P.60
% C
% C  CALCULATES COMPONENTS  OF REGION 1/2 FIELD IN SPHERICAL COORDS.  DERIVED FROM THE S/R DIPDEF2C (WHICH
% C    DOES THE SAME JOB, BUT INPUT/OUTPUT THERE WAS IN SPHERICAL COORDS, WHILE HERE WE USE CARTESIAN ONES)
% C
% C   INPUT:  NUMB=1 (2) FOR REGION 1 (2) CURRENTS
% C           MODE=1 YIELDS SIMPLE SINUSOIDAL MLT VARIATION, WITH MAXIMUM CURRENT AT DAWN/DUSK MERIDIAN
% C     WHILE MODE=2 YIELDS THE SECOND HARMONIC.
% C
% C
%      IMPLICIT REAL*8 (A-H,O-Z)
%      DIMENSION A11(31),A12(31),A21(31),A22(31)
%      COMMON /MODENUM/ M
%      COMMON /DTHETA/ DTHETA
%      COMMON /DPHI_B_RHO0/ DPHI,B,RHO_0,XKAPPA %! THESE PARAMETERS CONTROL DAY-NIGHT ASYMMETRY OF F.A.C., AS FOLLOWS:
global T01; % M and DTHETA are scalars in T01

% C  (1) DPHI:   HALF-DIFFERENCE (IN RADIANS) BETWEEN DAY AND NIGHT LATITUDE OF FAC OVAL AT IONOSPHERIC ALTITUDE;
% C              TYPICAL VALUE: 0.06
% C  (2) B:      AN ASYMMETRY FACTOR AT HIGH-ALTITUDES;  FOR B=0, THE ONLY ASYMMETRY IS THAT FROM DPHI
% C              TYPICAL VALUES: 0.35-0.70
% C  (3) RHO_0:  A FIXED PARAMETER, DEFINING THE DISTANCE RHO, AT WHICH THE LATITUDE SHIFT GRADUALLY SATURATES AND
% C              STOPS INCREASING
% C              ITS VALUE WAS ASSUMED FIXED, EQUAL TO 7.0.
% C  (4) XKAPPA: AN OVERALL SCALING FACTOR, WHICH CAN BE USED FOR CHANGING THE SIZE OF THE F.A.C. OVAL
% C
%      DATA BETA,RH,EPS/0.9D0,10.D0,3.D0/ %! parameters of the tilt-dependent deformation of the untilted F.A.C. field
BETA = 0.9D0;
RH = 10.D0;
EPS = 3.D0;

A11 = [.1618068350,-.1797957553,2.999642482,-.9322708978, ...  
        -.6811059760,.2099057262,-8.358815746,-14.86033550,.3838362986, ...
        -16.30945494,4.537022847,2.685836007,27.97833029,6.330871059, ...
        1.876532361,18.95619213,.9651528100,.4217195118,-.08957770020, ...
        -1.823555887,.7457045438,-.5785916524,-1.010200918,.01112389357, ...
        .09572927448,-.3599292276,8.713700514,.9763932955,3.834602998, ...
        2.492118385,.7113544659];
A12 = [.7058026940,-.2845938535,5.715471266,-2.472820880, ...
        -.7738802408,.3478293930,-11.37653694,-38.64768867,.6932927651, ...
        -212.4017288,4.944204937,3.071270411,33.05882281,7.387533799, ...
        2.366769108,79.22572682,.6154290178,.5592050551,-.1796585105, ...
        -1.654932210,.7309108776,-.4926292779,-1.130266095,-.009613974555, ...
        .1484586169,-.2215347198,7.883592948,.02768251655,2.950280953, ...
        1.212634762,.5567714182];
A21 = [.1278764024,-.2320034273,1.805623266,-32.37241440, ...
        -.9931490648,.3175085630,-2.492465814,-16.21600096,.2695393416, ...
        -6.752691265,3.971794901,14.54477563,41.10158386,7.912889730, ...
        1.258297372,9.583547721,1.014141963,.5104134759,-.1790430468, ...
        -1.756358428,.7561986717,-.6775248254,-.04014016420,.01446794851, ...
        .1200521731,-.2203584559,4.508963850,.8221623576,1.779933730, ...
        1.102649543,.8867880020];
A22 = [.4036015198,-.3302974212,2.827730930,-45.44405830, ...
        -1.611103927,.4927112073,-.003258457559,-49.59014949,.3796217108, ...
        -233.7884098,4.312666980,18.05051709,28.95320323,11.09948019, ...
        .7471649558,67.10246193,.5667096597,.6468519751,-.1560665317, ...
        -1.460805289,.7719653528,-.6658988668,.2515179349E-05, ...
        .02426021891,.1195003324,-.2625739255,4.377172556,.2421190547, ...
        2.503482679,1.071587299,.7247997430];

T01.DPHI_B_RHO0.B=0.5;
T01.DPHI_B_RHO0.RHO_0=7.0;

T01.M=MODE;
if (NUMB==1) ,
    T01.DPHI_B_RHO0.DPHI=0.055D0;
    T01.DTHETA=0.06D0;
end

if (NUMB==2) ,
    T01.DPHI_B_RHO0.DPHI=0.030D0;
    T01.DTHETA=0.09D0;
end

Xsc=X*T01.DPHI_B_RHO0.XKAPPA;
Ysc=Y*T01.DPHI_B_RHO0.XKAPPA;
Zsc=Z*T01.DPHI_B_RHO0.XKAPPA;
RHO=sqrt(Xsc^2+Zsc^2);

Rsc=sqrt(Xsc^2+Ysc^2+Zsc^2)  ;                               %!  SCALED
RHO2=T01.DPHI_B_RHO0.RHO_0^2;

if (Xsc==0.D0) & (Zsc==0.D0) ,
    PHI=0.D0;
else
    PHI=atan2(-Zsc,Xsc);  %!  FROM CARTESIAN TO CYLINDRICAL (RHO,PHI,Y)
end

SPHIC=sin(PHI);
CPHIC=cos(PHI);  %!  "C" means "CYLINDRICAL", TO DISTINGUISH FROM SPHERICAL PHI

BRACK=T01.DPHI_B_RHO0.DPHI+T01.DPHI_B_RHO0.B*RHO2/(RHO2+1.D0)*(RHO^2-1.D0)/(RHO2+RHO^2);
R1RH=(Rsc-1.D0)/RH;
if (R1RH<0.D0), R1RH=0.D0; end  %!  AVOID NEGATIVE VALUES OF R1RH, WHICH MAY OCCUR IF THE
% C                                     POINT (X,Y,Z) LIES CLOSE TO EARTH'S SURFACE AND THE S.W.
% C                                     PRESSURE IS ABNORMALLY LOW

PSIAS=BETA*PS/(1.D0+R1RH^EPS)^(1.D0/EPS);

PHIS=PHI-BRACK*sin(PHI) -PSIAS;
DPHISPHI=1.D0-BRACK*cos(PHI);
DPHISRHO=-2.D0*T01.DPHI_B_RHO0.B*RHO2*RHO/(RHO2+RHO^2)^2 *sin(PHI) ...
    +BETA*PS*R1RH^(EPS-1.D0)*RHO/(RH*Rsc* ...
    (1.D0+R1RH^EPS)^(1.D0/EPS+1.D0));
DPHISDY= BETA*PS*R1RH^(EPS-1.D0)*Ysc/(RH*Rsc* ...
    (1.D0+R1RH^EPS)^(1.D0/EPS+1.D0));

SPHICS=sin(PHIS);
CPHICS=cos(PHIS);

XS= RHO*CPHICS;
ZS=-RHO*SPHICS;

if (NUMB==1) ,
    if (MODE==1), [BXS,BYAS,BZS] = T01_TWOCONES(A11,XS,Ysc,ZS); end
    if (MODE==2), [BXS,BYAS,BZS] = T01_TWOCONES(A12,XS,Ysc,ZS); end
else
    if (MODE==1), [BXS,BYAS,BZS] = T01_TWOCONES(A21,XS,Ysc,ZS); end
    if (MODE==2), [BXS,BYAS,BZS] = T01_TWOCONES(A22,XS,Ysc,ZS); end
end

BRHOAS=BXS*CPHICS-BZS*SPHICS;
BPHIAS=-BXS*SPHICS-BZS*CPHICS;

BRHO_S=BRHOAS*DPHISPHI                             *T01.DPHI_B_RHO0.XKAPPA;        %! SCALING
BPHI_S=(BPHIAS-RHO*(BYAS*DPHISDY+BRHOAS*DPHISRHO)) *T01.DPHI_B_RHO0.XKAPPA;
BY_S=BYAS*DPHISPHI                                 *T01.DPHI_B_RHO0.XKAPPA;

BX=BRHO_S*CPHIC-BPHI_S*SPHIC;
BY=BY_S;
BZ=-BRHO_S*SPHIC-BPHI_S*CPHIC;

% end of function BIRK_1N2      

% c
% C=========================================================================
% c
function [BX,BY,BZ] = T01_TWOCONES (A,X,Y,Z);
%      SUBROUTINE TWOCONES (A,X,Y,Z,BX,BY,BZ)
% C
% C  ADDS FIELDS FROM TWO CONES (NORTHERN AND SOUTHERN), WITH A PROPER SYMMETRY
% C  OF THE CURRENT AND FIELD, CORRESPONDING TO THE REGION 1 BIRKELAND CURRENTS. (NB #6, P.58).
% C
%      IMPLICIT REAL*8 (A-H,O-Z)
%      DIMENSION A(31)

[BXN,BYN,BZN] = T01_ONE_CONE(A,X,Y,Z);
[BXS,BYS,BZS] = T01_ONE_CONE(A,X,-Y,-Z);
BX=BXN-BXS;
BY=BYN+BYS;
BZ=BZN+BZS;

% end of function TWOCONES

% c
% C-------------------------------------------------------------------------
% C
function [BX,BY,BZ] = T01_ONE_CONE(A,X,Y,Z);
%      SUBROUTINE ONE_CONE(A,X,Y,Z,BX,BY,BZ)
% c
% c  RETURNS FIELD COMPONENTS FOR A DEFORMED CONICAL CURRENT SYSTEM, FITTED TO A BIOSAVART FIELD
% c  HERE ONLY THE NORTHERN CONE IS TAKEN INTO ACCOUNT.
% c
%      IMPLICIT REAL*8 (A-H,O-Z)
%      DIMENSION A(31)

%      COMMON /DTHETA/ DTHETA
%      COMMON /MODENUM/ M
global T01;

%      DATA DR,DT/1.D-6,1.D-6/  %!   JUST FOR NUMERICAL DIFFERENTIATION
DR = 1.D-6;
DT =1.D-6;  %!   JUST FOR NUMERICAL DIFFERENTIATION

THETA0=A(31);

RHO2=X^2+Y^2;
RHO=sqrt(RHO2);
R=sqrt(RHO2+Z^2);
THETA=atan2(RHO,Z);
PHI=atan2(Y,X);
% C
% C   MAKE THE DEFORMATION OF COORDINATES:
% C
RS=T01_R_S(A,R,THETA);
THETAS=T01_THETA_S(A,R,THETA);
PHIS=PHI;
% C
% C   CALCULATE FIELD COMPONENTS AT THE NEW POSITION (ASTERISKED):
% C
[BTAST,BFAST] = T01_FIALCOS (RS,THETAS,PHIS,T01.M,THETA0,T01.DTHETA);    %!   MODE #M
% C
% C   NOW TRANSFORM B{R,T,F}_AST BY THE DEFORMATION TENSOR:
% C
% C      FIRST OF ALL, FIND THE DERIVATIVES:
% C
DRSDR=(T01_R_S(A,R+DR,THETA)-T01_R_S(A,R-DR,THETA))/(2.D0*DR);
DRSDT=(T01_R_S(A,R,THETA+DT)-T01_R_S(A,R,THETA-DT))/(2.D0*DT);
DTSDR=(T01_THETA_S(A,R+DR,THETA)-T01_THETA_S(A,R-DR,THETA))/(2.D0*DR);
DTSDT=(T01_THETA_S(A,R,THETA+DT)-T01_THETA_S(A,R,THETA-DT))/(2.D0*DT);

STSST=sin(THETAS)/sin(THETA);
RSR=RS/R;

BR     =-RSR/R*STSST*BTAST*DRSDT;                 %!   NB#6, P.43    BRAST DOES NOT ENTER HERE
BTHETA = RSR*STSST*BTAST*DRSDR;                  %!          (IT IS IDENTICALLY ZERO IN OUR CASE)
BPHI   = RSR*BFAST*(DRSDR*DTSDT-DRSDT*DTSDR);

S=RHO/R;
C=Z/R;
SF=Y/RHO;
CF=X/RHO;

BE=BR*S+BTHETA*C;

BX=A(1)*(BE*CF-BPHI*SF);
BY=A(1)*(BE*SF+BPHI*CF);
BZ=A(1)*(BR*C-BTHETA*S);

% end of function ONE_CONE

% C
% C=====================================================================================
function R_S = T01_R_S(A,R,THETA)
%      DOUBLE PRECISION FUNCTION R_S(A,R,THETA)
%      IMPLICIT REAL*8 (A-H,O-Z)
%      DIMENSION A(31)
% C
R_S=R+A(2)/R+A(3)*R/sqrt(R^2+A(11)^2)+A(4)*R/(R^2+A(12)^2) ...
    +(A(5)+A(6)/R+A(7)*R/sqrt(R^2+A(13)^2)+A(8)*R/(R^2+A(14)^2))* ...
    cos(THETA) ...
    +(A(9)*R/sqrt(R^2+A(15)^2)+A(10)*R/(R^2+A(16)^2)^2)* ...
    cos(2.D0*THETA);
    % C

% end of function R_S

% C
% C-----------------------------------------------------------------------------
% C
function THETA_S = T01_THETA_S(A,R,THETA)
%      DOUBLE PRECISION FUNCTION THETA_S(A,R,THETA)
%      IMPLICIT REAL*8 (A-H,O-Z)
%      DIMENSION A(31)
% c
THETA_S=THETA+(A(17)+A(18)/R+A(19)/R^2 ...
    +A(20)*R/sqrt(R^2+A(27)^2))*sin(THETA) ...
    +(A(21)+A(22)*R/sqrt(R^2+A(28)^2) ...
    +A(23)*R/(R^2+A(29)^2))*sin(2.D0*THETA) ...
    +(A(24)+A(25)/R+A(26)*R/(R^2+A(30)^2))*sin(3.D0*THETA);
    % C
% end of function THETA_S

% C
% c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
% c
function [BTHETA,BPHI] = T01_FIALCOS(R,THETA,PHI,N,THETA0,DT);
%      SUBROUTINE FIALCOS(R,THETA,PHI,BTHETA,BPHI,N,THETA0,DT)
% C
% C  CONICAL MODEL OF BIRKELAND CURRENT FIELD; BASED ON THE OLD S/R FIALCO (OF 1990-91)
% C  NB OF 1985-86-88, NOTE OF MARCH 5, BUT HERE BOTH INPUT AND OUTPUT ARE IN SPHERICAL CDS.
% 
% C  BTN, AND BPN ARE THE ARRAYS OF BTHETA AND BPHI (BTN(i), BPN(i) CORRESPOND TO i-th MODE).
% C   ONLY FIRST  N  MODE AMPLITUDES ARE COMPUTED (N&lt=10).
% C    THETA0 IS THE ANGULAR HALF-WIDTH OF THE CONE, DT IS THE ANGULAR H.-W. OF THE CURRENT LAYER
% 
% C   NOTE:  BR=0  (BECAUSE ONLY RADIAL CURRENTS ARE PRESENT IN THIS MODEL)
% C
%      IMPLICIT REAL*8 (A-H,O-Z)
%      DIMENSION  BTN(10),BPN(10),CCOS(10),SSIN(10)

SINTE=sin(THETA);
RO=R*SINTE;
COSTE=cos(THETA);
SINFI=sin(PHI);
COSFI=cos(PHI);
TG=SINTE/(1.D0+COSTE);   %!        TAN(THETA/2)
CTG=SINTE/(1.D0-COSTE);  %!        COT(THETA/2)
% C
% C
TETANP=THETA0+DT;
TETANM=THETA0-DT;
if (THETA>=TETANM),
    TGP=tan(TETANP*0.5D0);
    TGM=tan(TETANM*0.5D0);
    TGM2=TGM*TGM;
    TGP2=TGP*TGP;
end

COSM1=1.D0;
SINM1=0.D0;
TM=1.D0;
TGM2M=1.D0;
TGP2M=1.D0;

for M=1:N,
    TM=TM*TG;
    CCOS(M)=COSM1*COSFI-SINM1*SINFI;
    SSIN(M)=SINM1*COSFI+COSM1*SINFI;
    COSM1=CCOS(M);
    SINM1=SSIN(M);
    if (THETA<TETANM) ,
        T=TM;
        DTT=0.5D0*M*TM*(TG+CTG);
        DTT0=0.D0;
    elseif (THETA<TETANP) ,
        TGM2M=TGM2M*TGM2;
        FC=1.D0/(TGP-TGM);
        FC1=1.D0/(2*M+1);
        TGM2M1=TGM2M*TGM;
        TG21=1.D0+TG*TG;
        T=FC*(TM*(TGP-TG)+FC1*(TM*TG-TGM2M1/TM));
        DTT=0.5D0*M*FC*TG21*(TM/TG*(TGP-TG)-FC1*(TM-TGM2M1/(TM*TG)));
        DTT0=0.5D0*FC*((TGP+TGM)*(TM*TG-FC1*(TM*TG-TGM2M1/TM))+ ...
            TM*(1.D0-TGP*TGM)-(1.D0+TGM2)*TGM2M/TM);
    else
        TGP2M=TGP2M*TGP2;
        TGM2M=TGM2M*TGM2;
        FC=1.D0/(TGP-TGM);
        FC1=1.D0/(2*M+1);
        T=FC*FC1*(TGP2M*TGP-TGM2M*TGM)/TM;
        DTT=-T*M*0.5D0*(TG+CTG);
    end
    
    BTN(M)=M*T*CCOS(M)/RO;
    BPN(M)=-DTT*SSIN(M)/R;
end

BTHETA=BTN(N) *800.;
BPHI  =BPN(N) *800.;

% end of function FIALCOS

% C
% C-------------------------------------------------------------------------
% C
% C
function [BX,BY,BZ] = T01_BIRK_SHL(A,PS,X_SC,X,Y,Z);
%         SUBROUTINE BIRK_SHL (A,PS,X_SC,X,Y,Z,BX,BY,BZ)
% C
%	 IMPLICIT  REAL * 8  (A - H, O - Z)
%         DIMENSION A(86)
% C
CPS=cos(PS);
SPS=sin(PS);

S3PS=2.D0*CPS;
% C
PST1=PS*A(85);
PST2=PS*A(86);

ST1=sin(PST1);
CT1=cos(PST1);
ST2=sin(PST2);
CT2=cos(PST2);

X1=X*CT1-Z*ST1;
Z1=X*ST1+Z*CT1;
X2=X*CT2-Z*ST2;
Z2=X*ST2+Z*CT2;
% C
L=0;
GX=0.D0;
GY=0.D0;
GZ=0.D0;
% C
for M=1:2,
    %         DO 1 M=1,2     %!    M=1 IS FOR THE 1ST SUM ("PERP." SYMMETRY)
    % C                          AND M=2 IS FOR THE SECOND SUM ("PARALL." SYMMETRY)
    for I=1:3,
        %             DO 2 I=1,3
        P=A(72+I);
        Q=A(78+I);
        CYPI=cos(Y/P);
        CYQI=cos(Y/Q);
        SYPI=sin(Y/P);
        SYQI=sin(Y/Q);
        % C
        for K=1:3,
            %                DO 3 K=1,3
            R=A(75+K);
            S=A(81+K);
            SZRK=sin(Z1/R);
            CZSK=cos(Z2/S);
            CZRK=cos(Z1/R);
            SZSK=sin(Z2/S);
            SQPR=sqrt(1.D0/P^2+1.D0/R^2);
            SQQS=sqrt(1.D0/Q^2+1.D0/S^2);
            EPR=exp(X1*SQPR);
            EQS=exp(X2*SQQS);
            % C
            for N=1:2,
                %                  DO 4 N=1,2  %! N=1 IS FOR THE FIRST PART OF EACH COEFFICIENT
                % C                                AND N=2 IS FOR THE SECOND ONE
                for NN = 1:2,
                    %                    DO 5 NN=1,2 %!   NN = 1,2 FURTHER SPLITS THE COEFFICIENTS INTO 2 PARTS,
                    % C                                         TO TAKE INTO ACCOUNT THE SCALE FACTOR DEPENDENCE
                    
                    if (M==1),
                        FX=-SQPR*EPR*CYPI*SZRK;
                        FY=EPR*SYPI*SZRK/P;
                        FZ=-EPR*CYPI*CZRK/R;
                        if (N==1) ,
                            if (NN==1) ,
                                HX=FX;
                                HY=FY;
                                HZ=FZ;
                            else
                                HX=FX*X_SC;
                                HY=FY*X_SC;
                                HZ=FZ*X_SC;
                            end
                        else
                            if (NN==1) ,
                                HX=FX*CPS;
                                HY=FY*CPS;
                                HZ=FZ*CPS;
                            else
                                HX=FX*CPS*X_SC;
                                HY=FY*CPS*X_SC;
                                HZ=FZ*CPS*X_SC;
                            end
                        end
                        
                    else                            %!   M==2
                        FX=-SPS*SQQS*EQS*CYQI*CZSK;
                        FY=SPS/Q*EQS*SYQI*CZSK;
                        FZ=SPS/S*EQS*CYQI*SZSK;
                        if (N==1) ,
                            if (NN==1) ,
                                HX=FX;
                                HY=FY;
                                HZ=FZ;
                            else
                                HX=FX*X_SC;
                                HY=FY*X_SC;
                                HZ=FZ*X_SC;
                            end
                        else
                            if (NN==1) ,
                                HX=FX*S3PS;
                                HY=FY*S3PS;
                                HZ=FZ*S3PS;
                            else
                                HX=FX*S3PS*X_SC;
                                HY=FY*S3PS*X_SC;
                                HZ=FZ*S3PS*X_SC;
                            end
                        end
                    end
                    L=L+1;
                    
                    if (M==1) ,
                        HXR=HX*CT1+HZ*ST1;
                        HZR=-HX*ST1+HZ*CT1;
                    else
                        HXR=HX*CT2+HZ*ST2;
                        HZR=-HX*ST2+HZ*CT2;
                    end
                    
                    GX=GX+HXR*A(L);
                    GY=GY+HY *A(L);
                    GZ=GZ+HZR*A(L);
                end %5
                
            end %  4   CONTINUE
        end %  3   CONTINUE
    end %  2   CONTINUE
end %  1   CONTINUE

BX=GX;
BY=GY;
BZ=GZ;

% end of function BIRK_SHL

% C
% C^**********************************************************************************
% C
function [BXSRC,BYSRC,BZSRC,BXPRC,BYPRC,BZPRC] = T01_FULL_RC(IOPR,PS,X,Y,Z);
%      SUBROUTINE FULL_RC (IOPR,PS,X,Y,Z,BXSRC,BYSRC,BZSRC,BXPRC,BYPRC,
%     *  BZPRC)
% C
% C   CALCULATES GSM FIELD COMPONENTS OF THE SYMMETRIC (SRC) AND PARTIAL (PRC) COMPONENTS OF THE RING CURRENT
% C   SRC  PROVIDES A DEPRESSION OF -28 nT AT EARTH
% C   PRC  CORRESPONDS TO THE PRESSURE DIFFERENCE OF 2 nPa BETWEEN MIDNIGHT AND NOON RING CURRENT
% C             PARTICLE PRESSURE AND YIELDS A DEPRESSION OF -17 nT AT X=-6Re
% C
% C   SC_SY AND SC_PR ARE SCALING FACTORS FOR THE SYMMETRIC AND PARTIAL COMPONENTS:
% C          VALUES LARGER THAN 1 RESULT IN SPATIALLY LARGER CURRENTS
% C
% C   PHI IS THE ROTATION ANGLE IN RADIANS OF THE PARTIAL RING CURRENT (MEASURED FROM MIDNIGHT TOWARD DUSK)
% C
% C     IOPR -  A RING CURRENT CALCULATION FLAG (FOR LEAST-SQUARES FITTING ONLY):
% C             IOPR=0 - BOTH SRC AND PRC FIELDS ARE CALCULATED
% C             IOPR=1 - SRC ONLY
% C             IOPR=2 - PRC ONLY
% C
%        IMPLICIT REAL*8 (A-H,O-Z)
%        DIMENSION C_SY(86),C_PR(86)
%        COMMON /RCPAR/ SC_SY,SC_PR,PHI
%%%%%% original:COMMON /RCPAR/ SC_SY,SC_AS,PHI
% note SC_PR -< T01.RCPAR.SC_AS
global T01;
% C
C_SY = [1675.694858,1780.006388,-961.6082149,-1668.914259, ...
        -27.40437029,-107.4169670,27.76189943,92.89740503,-43.92949274, ...
        -403.6444072,6.167161865,298.2779761,-1680.779044,-1780.933039, ...
        964.1861088,1670.988659,27.48864650,107.7809519,-27.84600972, ...
        -93.20691865,44.28496784,404.4537249,-6.281958730,-298.6050952, ...
        -7.971914848,2.017383761,-1.492230168,-1.957411655,-.08525523181, ...
        -.3811813235,.08446716725,.3215044399,-.7141912767,-.9086294596, ...
        .2966677742,-.04736679933,-11.38731325,.1719795189,1.356233066, ...
        .8613438429,-.09143823092,-.2593979098,.04244838338,.06318383319, ...
        -.5861372726,-.03368780733,-.07104470269,-.06909052953, ...
        -60.18659631,-32.87563877,11.76450433,5.891673644,2.562360333, ...
        6.215377232,-1.273945165,-1.864704763,-5.394837143,-8.799382627, ...
        3.743066561,-.7649164511,57.09210569,32.61236511,-11.28688017, ...
        -5.849523392,-2.470635922,-5.961417272,1.230031099,1.793192595, ...
        5.383736074,8.369895153,-3.611544412,.7898988697,7.970609948, ...
        7.981216562,35.16822497,12.45651654,1.689755359,3.678712366, ...
        23.66117284,6.987136092,6.886678677,20.91245928,1.650064156, ...
        3.474068566,.3474715765,.6564043111];

C_PR = [-64820.58481,-63965.62048,66267.93413,135049.7504, ...
        -36.56316878,124.6614669,56.75637955,-87.56841077,5848.631425, ...
        4981.097722,-6233.712207,-10986.40188,68716.52057,65682.69473, ...
        -69673.32198,-138829.3568,43.45817708,-117.9565488,-62.14836263, ...
        79.83651604,-6211.451069,-5151.633113,6544.481271,11353.03491, ...
        23.72352603,-256.4846331,25.77629189,145.2377187,-4.472639098, ...
        -3.554312754,2.936973114,2.682302576,2.728979958,26.43396781, ...
        -9.312348296,-29.65427726,-247.5855336,-206.9111326,74.25277664, ...
        106.4069993,15.45391072,16.35943569,-5.965177750,-6.079451700, ...
        115.6748385,-35.27377307,-32.28763497,-32.53122151,93.74409310, ...
        84.25677504,-29.23010465,-43.79485175,-6.434679514,-6.620247951, ...
        2.443524317,2.266538956,-43.82903825,6.904117876,12.24289401, ...
        17.62014361,152.3078796,124.5505289,-44.58690290,-63.02382410, ...
        -8.999368955,-9.693774119,3.510930306,3.770949738,-77.96705716, ...
        22.07730961,20.46491655,18.67728847,9.451290614,9.313661792, ...
        644.7620970,418.2515954,7.183754387,35.62128817,19.43180682, ...
        39.57218411,15.69384715,7.123215241,2.300635346,21.90881131, ...
        -.01775839370,.3996346710];

[HXSRC,HYSRC,HZSRC,HXPRC,HYPRC,HZPRC] = T01_SRC_PRC(IOPR,...
    T01.RCPAR.SC_SY,T01.RCPAR.SC_AS,T01.RCPAR.PHI,PS,X,Y,Z);

X_SC=T01.RCPAR.SC_SY-1.D0;
if (IOPR==0) | (IOPR==1) ,
    [FSX,FSY,FSZ] = T01_RC_SHIELD (C_SY,PS,X_SC,X,Y,Z);
else
    FSX=0.D0;
    FSY=0.D0;
    FSZ=0.D0;
end

X_SC=T01.RCPAR.SC_AS-1.D0;
if (IOPR==0) | (IOPR==2) ,
    [FPX,FPY,FPZ] = T01_RC_SHIELD (C_PR,PS,X_SC,X,Y,Z);
else
    FPX=0.D0;
    FPY=0.D0;
    FPZ=0.D0;
end

BXSRC=HXSRC+FSX;
BYSRC=HYSRC+FSY;
BZSRC=HZSRC+FSZ;

BXPRC=HXPRC+FPX;
BYPRC=HYPRC+FPY;
BZPRC=HZPRC+FPZ;


% end of function FULL_RC

% C---------------------------------------------------------------------------------------
% C
function [BXSRC,BYSRC,BZSRC,BXPRC,BYPRC,BZPRC] = T01_SRC_PRC (IOPR,SC_SY,SC_PR,PHI,PS,X,Y,Z);
%        SUBROUTINE SRC_PRC (IOPR,SC_SY,SC_PR,PHI,PS,X,Y,Z,BXSRC,BYSRC,
%      *    BZSRC,BXPRC,BYPRC,BZPRC)
% C
% C   RETURNS FIELD COMPONENTS FROM A MODEL RING CURRENT, INCLUDING ITS SYMMETRIC PART
% C     AND A PARTIAL RING CURRENT, CLOSED VIA BIRKELAND CURRENTS. BASED ON RESULTS, DESCRIBED
% C     IN A PAPER "MODELING THE INNER MAGNETOSPHERE: ASYMMETRIC RING CURRENT AND REGION 2
% C     BIRKELAND CURRENTS REVISITED" (JGR, DEC.2000).
% C
% C     IOPR -  A RING CURRENT CALCULATION FLAG (FOR LEAST-SQUARES FITTING ONLY):
% C             IOPR=0 - BOTH SRC AND PRC FIELDS ARE CALCULATED
% C             IOPR=1 - SRC ONLY
% C             IOPR=2 - PRC ONLY
% C
% C     SC_SY &  SC_PR ARE SCALE FACTORS FOR THE ABOVE COMPONENTS;  TAKING SC<1 OR SC>1 MAKES THE CURRENTS
% C                      SHRINK OR EXPAND, RESPECTIVELY.
% C
% C   PHI IS THE ROTATION ANGLE (RADIANS) OF THE PARTIAL RING CURRENT (MEASURED FROM MIDNIGHT TOWARD DUSK)
% C
%        IMPLICIT REAL*8 (A-H,O-Z)
% c
% c   1.  TRANSFORM TO TILTED COORDINATES (i.e., SM coordinates):
% C
CPS=cos(PS);
SPS=sin(PS);

XT=X*CPS-Z*SPS;
ZT=Z*CPS+X*SPS;
% C
% C   2.  SCALE THE COORDINATES FOR THE SYMMETRIC AND PARTIAL RC COMPONENTS:
% C
XTS=XT/SC_SY;    %!  SYMMETRIC
YTS=Y /SC_SY;
ZTS=ZT/SC_SY;

XTA=XT/SC_PR;    %!  PARTIAL
YTA=Y /SC_PR;
ZTA=ZT/SC_PR;
% C
% C   3.  CALCULATE COMPONENTS OF THE TOTAL FIELD IN THE TILTED (SOLAR-MAGNETIC) COORDINATE SYSTEM:
% C
% C==========   ONLY FOR LEAST SQUARES FITTING:
BXS=0.D0;
BYS=0.D0;
BZS=0.D0;
BXA_S=0.D0;
BYA_S=0.D0;
BZA_S=0.D0;
BXA_QR=0.D0;
BYA_QR=0.D0;
BZA_Q=0.D0;
% C============================================
% C
% C    3a. SYMMETRIC FIELD:
% C
if (IOPR<=1), [BXS,BYS,BZS] = T01_RC_SYMM(XTS,YTS,ZTS); end;
if (IOPR==0) | (IOPR==2),
    [BXA_S,BYA_S,BZA_S] = T01_PRC_SYMM(XTA,YTA,ZTA);
end

% C    3b. ROTATE THE SCALED SM COORDINATES BY PHI AROUND ZSM AXIS AND CALCULATE QUADRUPOLE PRC FIELD
% C         IN THOSE COORDS:

CP=cos(PHI);
SP=sin(PHI);
XR=XTA*CP-YTA*SP;
YR=XTA*SP+YTA*CP;

if (IOPR==0) | (IOPR==2),
    [BXA_QR,BYA_QR,BZA_Q] = T01_PRC_QUAD(XR,YR,ZTA);
end

% C    3c. TRANSFORM THE QUADRUPOLE FIELD COMPONENTS BACK TO THE SM COORDS:
% C
BXA_Q= BXA_QR*CP+BYA_QR*SP;
BYA_Q=-BXA_QR*SP+BYA_QR*CP;

% C    3d. FIND THE TOTAL FIELD OF PRC (SYMM.+QUADR.) IN THE SM COORDS:
% C
BXP=BXA_S+BXA_Q;
BYP=BYA_S+BYA_Q;
BZP=BZA_S+BZA_Q;
% C
% C   4.  TRANSFORM THE FIELDS OF BOTH PARTS OF THE RING CURRENT BACK TO THE GSM SYSTEM:
% C
BXSRC=BXS*CPS+BZS*SPS;   %!    SYMMETRIC RC
BYSRC=BYS;
BZSRC=BZS*CPS-BXS*SPS;
% C
BXPRC=BXP*CPS+BZP*SPS;   %!    PARTIAL RC
BYPRC=BYP;
BZPRC=BZP*CPS-BXP*SPS;
% C

% end of function SRC_PRC

% C
% C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
% C
function [BX,BY,BZ]= T01_RC_SYMM(X,Y,Z)
%      SUBROUTINE RC_SYMM (X,Y,Z,BX,BY,BZ)
%      IMPLICIT  REAL * 8  (A - H, O - Z)
%      DATA DS,DC/1.D-2,0.99994999875D0/, D/1.D-4/,DRD/5.D3/  %! DS=SIN(THETA) AT THE BOUNDARY OF THE LINEARITY
%          REGION; DC=sqrt(1-DS^2);  DRD=1/(2*D)
%! DS=SIN(THETA) AT THE BOUNDARY OF THE LINEARITY
DS = 1.D-2;
DC = 0.99994999875D0;
D = 1.D-4;
DRD = 5.D3;
RHO2=X^2+Y^2;
R2=RHO2+Z^2;
R=sqrt(R2);
RP=R+D;
RM=R-D;
SINT=sqrt(RHO2)/R;
COST=Z/R;

if (SINT<DS) ,  %!  TOO CLOSE TO THE Z-AXIS; USING A LINEAR APPROXIMATION A_PHI~SINT,
    % C                                    TO AVOID THE SINGULARITY PROBLEM
    A=T01_AP(R,DS,DC)/DS;
    DARDR=(RP*T01_AP(RP,DS,DC)-RM*T01_AP(RM,DS,DC))*DRD;
    FXY=Z*(2.D0*A-DARDR)/(R*R2);
    BX=FXY*X;
    BY=FXY*Y;
    BZ=(2.D0*A*COST^2+DARDR*SINT^2)/R;
    
else
    
    THETA=atan2(SINT,COST);
    TP=THETA+D;
    TM=THETA-D;
    SINTP=sin(TP);
    SINTM=sin(TM);
    COSTP=cos(TP);
    COSTM=cos(TM);
    BR=(SINTP*T01_AP(R,SINTP,COSTP)-SINTM*T01_AP(R,SINTM,COSTM)) ...
        /(R*SINT)*DRD;
    BT=(RM*T01_AP(RM,SINT,COST)-RP*T01_AP(RP,SINT,COST))/R*DRD;
    FXY=(BR+BT*COST/SINT)/R;
    BX=FXY*X;
    BY=FXY*Y;
    BZ=BR*COST-BT*SINT;
    
end

% end of function RC_SYMM

% c
% c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
% C
function AP = T01_AP(R,SINT,COST)
%      DOUBLE PRECISION FUNCTION AP(R,SINT,COST)
% C
% C      Calculates azimuthal component of the vector potential of the symmetric
% c  part of the model ring current.
% C
%      IMPLICIT  REAL * 8  (A - H, O - Z)
%      LOGICAL PROX   %!  INDICATES WHETHER WE ARE TOO CLOSE TO THE AXIS OF SYMMETRY, WHERE THE INVERSION
% C                                                             OF DIPOLAR COORDINATES BECOMES INACCURATE
%       DATA A1,A2,RRC1,DD1,RRC2,DD2,P1,R1,DR1,DLA1,P2,R2,DR2,DLA2,P3,
%      *R3,DR3/-563.3722359,425.0891691,4.150588549,2.266150226,
%      * 3.334503403,3.079071195,.02602428295,8.937790598,3.327934895,
%      *.4487061833,.09125832351,6.243029867,1.750145910,.4181957162,
%      *.06106691992,2.079908581,.6828548533/
% 
A1 = -563.3722359;
A2 = 425.0891691;
RRC1 = 4.150588549;
DD1 = 2.266150226;
RRC2 = 3.334503403;
DD2 = 3.079071195;
P1 = .02602428295;
R1 = 8.937790598;
DR1 = 3.327934895;
DLA1 = .4487061833;
P2 = .09125832351;
R2 = 6.243029867;
DR2 = 1.750145910;
DLA2 = .4181957162;
P3 = .06106691992;
R3 = 2.079908581;
DR3 = .6828548533;     

PROX=0;
SINT1=SINT;
COST1=COST;
if (SINT1<1.D-2) ,  %!  TOO CLOSE TO Z-AXIS;  USE LINEAR INTERPOLATION BETWEEN SINT=0 & SINT=0.01
    SINT1=1.D-2;
    COST1=.99994999875;
    PROX=1;
end

ALPHA=SINT1^2/R;         %!  R,THETA -> ALPHA,GAMMA
GAMMA=COST1/R^2;

ARG1=-((R-R1)/DR1)^2-(COST1/DLA1)^2;
ARG2=-((R-R2)/DR2)^2-(COST1/DLA2)^2;
ARG3=-((R-R3)/DR3)^2;

if (ARG1<-500.D0) ,        %!   TO PREVENT "FLOATING UNDERFLOW" CRASHES
    DEXP1=0.D0;
else
    DEXP1=exp(ARG1);
end

if (ARG2<-500.D0) ,
    DEXP2=0.D0;
else
    DEXP2=exp(ARG2);
end

if (ARG3<-500.D0) ,
    DEXP3=0.D0;
else
    DEXP3=exp(ARG3);
end


ALPHA_S=ALPHA*(1.D0+P1*DEXP1+P2*DEXP2+P3*DEXP3);     %!  ALPHA -> ALPHA_S  (DEFORMED)

GAMMA_S=GAMMA;
GAMMAS2=GAMMA_S^2;


ALSQH=ALPHA_S^2/2.D0;            %!  ALPHA_S,GAMMA_S -> RS,SINTS,COSTS
F=64.D0/27.D0*GAMMAS2+ALSQH^2;
Q=(sqrt(F)+ALSQH)^(1.D0/3.D0);
C=Q-4.D0*GAMMAS2^(1.D0/3.D0)/(3.D0*Q);
if (C<0.D0), C=0.D0; end
G=sqrt(C^2+4.D0*GAMMAS2^(1.D0/3.D0));
RS=4.D0/((sqrt(2.D0*G-C)+sqrt(C))*(G+C));
COSTS=GAMMA_S*RS^2;
SINTS=sqrt(1.D0-COSTS^2);
RHOS=RS*SINTS;
RHOS2=RHOS^2;
ZS=RS*COSTS;
% C
% c  1st loop:

P=(RRC1+RHOS)^2+ZS^2+DD1^2;
XK2=4.D0*RRC1*RHOS/P;
XK=sqrt(XK2);
XKRHO12=XK*sqrt(RHOS);     %!   SEE NB#4, P.3
% C
XK2S=1.D0-XK2;
DL=log(1.D0/XK2S);
ELK=1.38629436112d0+XK2S*(0.09666344259D0+XK2S*(0.03590092383+ ...
    XK2S*(0.03742563713+XK2S*0.01451196212))) +DL* ...
    (0.5D0+XK2S*(0.12498593597D0+XK2S*(0.06880248576D0+ ...
    XK2S*(0.03328355346D0+XK2S*0.00441787012D0))));
ELE=1.D0+XK2S*(0.44325141463D0+XK2S*(0.0626060122D0+XK2S* ...
    (0.04757383546D0+XK2S*0.01736506451D0))) +DL* ...
    XK2S*(0.2499836831D0+XK2S*(0.09200180037D0+XK2S* ...
    (0.04069697526D0+XK2S*0.00526449639D0)));
% C
APHI1=((1.D0-XK2*0.5D0)*ELK-ELE)/XKRHO12;
% c
% c  2nd loop:

P=(RRC2+RHOS)^2+ZS^2+DD2^2;
XK2=4.D0*RRC2*RHOS/P;
XK=sqrt(XK2);
XKRHO12=XK*sqrt(RHOS);     %!   SEE NB#4, P.3
% C
XK2S=1.D0-XK2;
DL=log(1.D0/XK2S);
ELK=1.38629436112d0+XK2S*(0.09666344259D0+XK2S*(0.03590092383+ ...
    XK2S*(0.03742563713+XK2S*0.01451196212))) +DL* ...
    (0.5D0+XK2S*(0.12498593597D0+XK2S*(0.06880248576D0+ ...
    XK2S*(0.03328355346D0+XK2S*0.00441787012D0))));
ELE=1.D0+XK2S*(0.44325141463D0+XK2S*(0.0626060122D0+XK2S* ...
    (0.04757383546D0+XK2S*0.01736506451D0))) +DL* ...
    XK2S*(0.2499836831D0+XK2S*(0.09200180037D0+XK2S* ...
    (0.04069697526D0+XK2S*0.00526449639D0)));
% C
APHI2=((1.D0-XK2*0.5D0)*ELK-ELE)/XKRHO12;

AP=A1*APHI1+A2*APHI2;
if (PROX), AP=AP*SINT/SINT1; end   %!   LINEAR INTERPOLATION, IF TOO CLOSE TO THE Z-AXIS
% C

% end of function AP


% c
% c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
% C
function [BX,BY,BZ] = T01_PRC_SYMM(X,Y,Z);
%      SUBROUTINE PRC_SYMM (X,Y,Z,BX,BY,BZ)
%      IMPLICIT  REAL * 8  (A - H, O - Z)
%! DS=SIN(THETA) AT THE BOUNDARY OF THE LINEARITY
%      REGION; DC=sqrt(1-DS^2);  DRD=1/(2*D)
DS = 1.D-2;
DC = 0.99994999875D0;
D = 1.D-4;
DRD = 5.D3;
RHO2=X^2+Y^2;
R2=RHO2+Z^2;
R=sqrt(R2);
RP=R+D;
RM=R-D;
SINT=sqrt(RHO2)/R;
COST=Z/R;

if (SINT<DS) ,  %!  TOO CLOSE TO THE Z-AXIS; USING A LINEAR APPROXIMATION A_PHI~SINT,
    % C                                    TO AVOID THE SINGULARITY PROBLEM
    A=T01_APPRC(R,DS,DC)/DS;
    DARDR=(RP*T01_APPRC(RP,DS,DC)-RM*T01_APPRC(RM,DS,DC))*DRD;
    FXY=Z*(2.D0*A-DARDR)/(R*R2);
    BX=FXY*X;
    BY=FXY*Y;
    BZ=(2.D0*A*COST^2+DARDR*SINT^2)/R;
    
else
    
    THETA=atan2(SINT,COST);
    TP=THETA+D;
    TM=THETA-D;
    SINTP=sin(TP);
    SINTM=sin(TM);
    COSTP=cos(TP);
    COSTM=cos(TM);
    BR=(SINTP*T01_APPRC(R,SINTP,COSTP)-SINTM*T01_APPRC(R,SINTM,COSTM)) ...
        /(R*SINT)*DRD;
    BT=(RM*T01_APPRC(RM,SINT,COST)-RP*T01_APPRC(RP,SINT,COST))/R*DRD;
    FXY=(BR+BT*COST/SINT)/R;
    BX=FXY*X;
    BY=FXY*Y;
    BZ=BR*COST-BT*SINT;
    
end

% end of function PRC_SYMM

% c
% c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
% C
% C
function APPRC = T01_APPRC(R,SINT,COST)
%      DOUBLE PRECISION FUNCTION APPRC(R,SINT,COST)
% C
% C      Calculates azimuthal component of the vector potential of the symmetric
% c  part of the model PARTIAL ring current.
% C
%      IMPLICIT  REAL * 8  (A - H, O - Z)
%      LOGICAL PROX
%       DATA A1,A2,RRC1,DD1,RRC2,DD2,P1,ALPHA1,DAL1,BETA1,DG1,P2,ALPHA2,
%      * DAL2,BETA2,DG2,BETA3,P3,ALPHA3,DAL3,BETA4,DG3,BETA5,Q0,Q1,ALPHA4,
%      * DAL4,DG4,Q2,ALPHA5,DAL5,DG5,BETA6,BETA7
%      * /-80.11202281,12.58246758,6.560486035,1.930711037,3.827208119,
%      *.7789990504,.3058309043,.1817139853,.1257532909,3.422509402,
%      *.04742939676,-4.800458958,-.02845643596,.2188114228,2.545944574,
%      *.00813272793,.35868244,103.1601001,-.00764731187,.1046487459,
%      *2.958863546,.01172314188,.4382872938,.01134908150,14.51339943,
%      *.2647095287,.07091230197,.01512963586,6.861329631,.1677400816,
%      *.04433648846,.05553741389,.7665599464,.7277854652/

A1 = -80.11202281;
A2 = 12.58246758;
RRC1 = 6.560486035;
DD1 = 1.930711037;
RRC2 = 3.827208119;
DD2 = .7789990504;
P1 = .3058309043;
ALPHA1 = .1817139853;
DAL1 = .1257532909;
BETA1 = 3.422509402;
DG1 = .04742939676;
P2 = -4.800458958;
ALPHA2 = -.02845643596;
DAL2 = .2188114228;
BETA2 = 2.545944574;
DG2 = .00813272793;
BETA3 = .35868244;
P3 = 103.1601001;
ALPHA3 = -.00764731187;
DAL3 = .1046487459;
BETA4 = 2.958863546;
DG3 = .01172314188;
BETA5 = .4382872938;
Q0 = .01134908150;
Q1 = 14.51339943;
ALPHA4 = .2647095287;
DAL4 = .07091230197;
DG4 = .01512963586;
Q2 = 6.861329631;
ALPHA5 = .1677400816;
DAL5 = .04433648846;
DG5 = .05553741389;
BETA6 = .7665599464;
BETA7 = .7277854652;

PROX=0;
SINT1=SINT;
COST1=COST;
if (SINT1<1.D-2) ,  %!  TOO CLOSE TO Z-AXIS;  USE LINEAR INTERPOLATION BETWEEN SINT=0 & SINT=0.01
    SINT1=1.D-2;
    COST1=.99994999875;
    PROX=1
end

ALPHA=SINT1^2/R;         %!  R,THETA -> ALPHA,GAMMA
GAMMA=COST1/R^2;

ALPHA_S=ALPHA*(1.D0+P1/(1.D0+((ALPHA-ALPHA1)/DAL1)^2)^BETA1 ...
    *exp(-(GAMMA/DG1)^2) ...
    +P2*(ALPHA-ALPHA2)/(1.D0+((ALPHA-ALPHA2)/DAL2)^2)^BETA2 ...
    /(1.D0+(GAMMA/DG2)^2)^BETA3 ...
    +P3*(ALPHA-ALPHA3)^2/(1.D0+((ALPHA-ALPHA3)/DAL3)^2)^BETA4 ...
    /(1.D0+(GAMMA/DG3)^2)^BETA5);     %!  ALPHA -> ALPHA_S  (DEFORMED)

GAMMA_S=GAMMA*(1.D0+Q0+Q1*(ALPHA-ALPHA4) ...
    *exp(-((ALPHA-ALPHA4)/DAL4)^2-(GAMMA/DG4)^2) ...   %!  GAMMA -> GAMMA_  (DEFORMED)
    +Q2*(ALPHA-ALPHA5)/(1.D0+((ALPHA-ALPHA5)/DAL5)^2)^BETA6 ...
    /(1.D0+(GAMMA/DG5)^2)^BETA7);

GAMMAS2=GAMMA_S^2;

ALSQH=ALPHA_S^2/2.D0;                            %!  ALPHA_S,GAMMA_S -> RS,SINTS,COSTS
F=64.D0/27.D0*GAMMAS2+ALSQH^2;
Q=(sqrt(F)+ALSQH)^(1.D0/3.D0);
C=Q-4.D0*GAMMAS2^(1.D0/3.D0)/(3.D0*Q);
if (C<0.D0), C=0.D0; end
G=sqrt(C^2+4.D0*GAMMAS2^(1.D0/3.D0));
RS=4.D0/((sqrt(2.D0*G-C)+sqrt(C))*(G+C));
COSTS=GAMMA_S*RS^2;
SINTS=sqrt(1.D0-COSTS^2);
RHOS=RS*SINTS;
RHOS2=RHOS^2;
ZS=RS*COSTS;
% C
% c  1st loop:

P=(RRC1+RHOS)^2+ZS^2+DD1^2;
XK2=4.D0*RRC1*RHOS/P;
XK=sqrt(XK2);
XKRHO12=XK*sqrt(RHOS);     %!    NB#4, P.3
% C
XK2S=1.D0-XK2;
DL=log(1.D0/XK2S);
ELK=1.38629436112d0+XK2S*(0.09666344259D0+XK2S*(0.03590092383+ ...
    XK2S*(0.03742563713+XK2S*0.01451196212))) +DL* ...
    (0.5D0+XK2S*(0.12498593597D0+XK2S*(0.06880248576D0+ ...
    XK2S*(0.03328355346D0+XK2S*0.00441787012D0))));
ELE=1.D0+XK2S*(0.44325141463D0+XK2S*(0.0626060122D0+XK2S* ...
    (0.04757383546D0+XK2S*0.01736506451D0))) +DL* ...
    XK2S*(0.2499836831D0+XK2S*(0.09200180037D0+XK2S* ...
    (0.04069697526D0+XK2S*0.00526449639D0)));
% C
APHI1=((1.D0-XK2*0.5D0)*ELK-ELE)/XKRHO12;
% c
% c  2nd loop:

P=(RRC2+RHOS)^2+ZS^2+DD2^2;
XK2=4.D0*RRC2*RHOS/P;
XK=sqrt(XK2);
XKRHO12=XK*sqrt(RHOS);     %!    NB#4, P.3
% C
XK2S=1.D0-XK2;
DL=log(1.D0/XK2S);
ELK=1.38629436112d0+XK2S*(0.09666344259D0+XK2S*(0.03590092383+ ...
    XK2S*(0.03742563713+XK2S*0.01451196212))) +DL* ...
    (0.5D0+XK2S*(0.12498593597D0+XK2S*(0.06880248576D0+ ...
    XK2S*(0.03328355346D0+XK2S*0.00441787012D0))));
ELE=1.D0+XK2S*(0.44325141463D0+XK2S*(0.0626060122D0+XK2S* ...
    (0.04757383546D0+XK2S*0.01736506451D0))) +DL* ...
    XK2S*(0.2499836831D0+XK2S*(0.09200180037D0+XK2S* ...
    (0.04069697526D0+XK2S*0.00526449639D0)));
% C
APHI2=((1.D0-XK2*0.5D0)*ELK-ELE)/XKRHO12;

APPRC=A1*APHI1+A2*APHI2;
if (PROX), APPRC=APPRC*SINT/SINT1; end   %!   LINEAR INTERPOLATION, IF TOO CLOSE TO THE Z-AXIS
% C
% end of function APPRC

% C
% C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
% C
% C
function [BX,BY,BZ] = T01_PRC_QUAD(X,Y,Z);
%         SUBROUTINE PRC_QUAD (X,Y,Z,BX,BY,BZ)
% C
% C  CALCULATES COMPONENTS OF THE FIELD FROM THE "QUADRUPOLE" COMPONENT OF THE PRC
% C
%         IMPLICIT  REAL * 8  (A - H, O - Z)

%         DATA D,DD/1.D-4,2.D-4/, DS/1.D-2/,DC/0.99994999875D0/
D = 1.D-4;
DD =2.D-4;
DS = 1.D-2;
DC = 0.99994999875D0;

RHO2=X^2+Y^2;
R=sqrt(RHO2+Z^2);
RHO=sqrt(RHO2);
SINT=RHO/R;
COST=Z/R;
RP=R+D;
RM=R-D;

if (SINT>DS) ,
    CPHI=X/RHO;
    SPHI=Y/RHO;
    BR=T01_BR_PRC_Q(R,SINT,COST);
    BT=T01_BT_PRC_Q(R,SINT,COST);
    DBRR=(T01_BR_PRC_Q(RP,SINT,COST)-T01_BR_PRC_Q(RM,SINT,COST))/DD;
    THETA=atan2(SINT,COST);
    TP=THETA+D;
    TM=THETA-D;
    SINTP=sin(TP);
    COSTP=cos(TP);
    SINTM=sin(TM);
    COSTM=cos(TM);
    DBTT=(T01_BT_PRC_Q(R,SINTP,COSTP)-T01_BT_PRC_Q(R,SINTM,COSTM))/DD;
    BX=SINT*(BR+(BR+R*DBRR+DBTT)*SPHI^2)+COST*BT;
    BY=-SINT*SPHI*CPHI*(BR+R*DBRR+DBTT);
    BZ=(BR*COST-BT*SINT)*CPHI;
else
    ST=DS;
    CT=DC;
    if (Z<0.D0), CT=-DC; end;
    THETA=atan2(ST,CT);
    TP=THETA+D;
    TM=THETA-D;
    SINTP=sin(TP);
    COSTP=cos(TP);
    SINTM=sin(TM);
    COSTM=cos(TM);
    BR=T01_BR_PRC_Q(R,ST,CT);
    BT=T01_BT_PRC_Q(R,ST,CT);
    DBRR=(T01_BR_PRC_Q(RP,ST,CT)-T01_BR_PRC_Q(RM,ST,CT))/DD;
    DBTT=(T01_BT_PRC_Q(R,SINTP,COSTP)-T01_BT_PRC_Q(R,SINTM,COSTM))/DD;
    FCXY=R*DBRR+DBTT;
    BX=(BR*(X^2+2.D0*Y^2)+FCXY*Y^2)/(R*ST)^2+BT*COST;
    BY=-(BR+FCXY)*X*Y/(R*ST)^2;
    BZ=(BR*COST/ST-BT)*X/R;
end
% end of function PRC_QUAD

% c
% c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
% C
function BR_PRC_Q = T01_BR_PRC_Q(R,SINT,COST)
%      DOUBLE PRECISION FUNCTION BR_PRC_Q (R,SINT,COST)
% C
% Calculates the radial component of the "quadrupole" part of the model partial ring current.
% C
%      IMPLICIT  REAL * 8  (A - H, O - Z)

%       DATA A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11,A12,A13,A14,A15,A16,A17,   %! ALL LINEAR PARAMETERS HERE
%      * A18,XK1,AL1,DAL1,B1,BE1,XK2,AL2,DAL2,B2,BE2,XK3,XK4,AL3,DAL3,B3,  %! WERE MULTIPLIED BY 0.1,
%      * BE3,AL4,DAL4,DG1,AL5,DAL5,DG2,C1,C2,C3,AL6,DAL6,DRM/-21.2666329,  %! SO THAT THEY CORRESPOND TO P_0=1 nPa,
%      *32.24527521,-6.062894078,7.515660734,233.7341288,-227.1195714,     %! RATHER THAN THE ORIGINAL VALUE OF 10 nPa
%      *8.483233889,16.80642754,-24.63534184,9.067120578,-1.052686913,     %! ASSUMED IN THE BIOT-SAVART INTEGRAL
%      *-12.08384538,18.61969572,-12.71686069,47017.35679,-50646.71204,
%      *7746.058231,1.531069371,2.318824273,.1417519429,.6388013110E-02,
%      *5.303934488,4.213397467,.7955534018,.1401142771,.2306094179E-01,
%      *3.462235072,2.568743010,3.477425908,1.922155110,.1485233485,
%      *.2319676273E-01,7.830223587,8.492933868,.1295221828,.01753008801,
%      *.01125504083,.1811846095,.04841237481,.01981805097,6.557801891,
%      *6.348576071,5.744436687,.2265212965,.1301957209,.5654023158/

A1 = -21.2666329;
A2 = 32.24527521;
A3 = -6.062894078;
A4 = 7.515660734;
A5 = 233.7341288;
A6 = -227.1195714;
A7 = 8.483233889;
A8 = 16.80642754;
A9 = -24.63534184;
A10 = 9.067120578;
A11 = -1.052686913;
A12 = -12.08384538;
A13 = 18.61969572;
A14 = -12.71686069;
A15 = 47017.35679;
A16 = -50646.71204;
A17 = 7746.058231;
A18 = 1.531069371;
XK1 = 2.318824273;
AL1 = .1417519429;
DAL1 = .6388013110E-02;
B1 = 5.303934488;
BE1 = 4.213397467;
XK2 = .7955534018;
AL2 = .1401142771;
DAL2 = .2306094179E-01;
B2 = 3.462235072;
BE2 = 2.568743010;
XK3 = 3.477425908;
XK4 = 1.922155110;
AL3 = .1485233485;
DAL3 = .2319676273E-01;
B3 = 7.830223587;
BE3 = 8.492933868;
AL4 = .1295221828;
DAL4 = .01753008801;
DG1 = .01125504083;
AL5 = .1811846095;
DAL5 = .04841237481;
DG2 = .01981805097;
C1 = 6.557801891;
C2 = 6.348576071;
C3 = 5.744436687;
AL6 = .2265212965;
DAL6 = .1301957209;
DRM = .5654023158;

SINT2=SINT^2;
COST2=COST^2;
SC=SINT*COST;
ALPHA=SINT2/R;
GAMMA=COST/R^2;
[F,FA,FS] = T01_FFS(ALPHA,AL1,DAL1);      
D1=SC*F^XK1/((R/B1)^BE1+1.D0);
D2=D1*COST2;

[F,FA,FS] = T01_FFS(ALPHA,AL2,DAL2);
D3=SC*FS^XK2/((R/B2)^BE2+1.D0);
D4=D3*COST2;

[F,FA,FS] = T01_FFS(ALPHA,AL3,DAL3);
D5=SC*(ALPHA^XK3)*(FS^XK4)/((R/B3)^BE3+1.D0);
D6=D5*COST2;

ARGA=((ALPHA-AL4)/DAL4)^2+1.D0;
ARGG=1.D0+(GAMMA/DG1)^2;

D7=SC/ARGA/ARGG;
D8=D7/ARGA;
D9=D8/ARGA;
D10=D9/ARGA;

ARGA=((ALPHA-AL5)/DAL5)^2+1.D0;
ARGG=1.D0+(GAMMA/DG2)^2;

D11=SC/ARGA/ARGG;
D12=D11/ARGA;
D13=D12/ARGA;
D14=D13/ARGA;


D15=SC/(R^4+C1^4);
D16=SC/(R^4+C2^4)*COST2;
D17=SC/(R^4+C3^4)*COST2^2;

[F,FA,FS] = T01_FFS(ALPHA,AL6,DAL6);
D18=SC*FS/(1.D0+((R-1.2D0)/DRM)^2);

BR_PRC_Q=A1*D1+A2*D2+A3*D3+A4*D4+A5*D5+A6*D6+A7*D7+A8*D8+A9*D9+ ...
    A10*D10+A11*D11+A12*D12+A13*D13+A14*D14+A15*D15+A16*D16+A17*D17+ ...
    A18*D18;
% C
% end of function BR_PRC_Q

% c
% C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% C
function BT_PRC_Q = T01_BT_PRC_Q(R,SINT,COST);
%        DOUBLE PRECISION FUNCTION BT_PRC_Q (R,SINT,COST)
% C
% Calculates the Theta component of the "quadrupole" part of the model partial ring current.
% C
%        IMPLICIT  REAL * 8  (A - H, O - Z)

%       DATA A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11,A12,A13,A14,A15,A16,A17,  %! ALL LINEAR PARAMETERS HERE
%      *XK1,AL1,DAL1,B1,BE1,XK2,AL2,DAL2,BE2,XK3,XK4,AL3,DAL3,B3,BE3,AL4, %! WERE MULTIPLIED BY 0.1,
%      *DAL4,DG1,AL5,DAL5,DG2,C1,C2,C3/12.74640393,-7.516393516,          %! SO THAT THEY CORRESPOND TO P_0=1 nPa,
%      *-5.476233865,3.212704645,-59.10926169,46.62198189,-.01644280062,  %! RATHER THAN THE ORIGINAL VALUE OF 10 nPa
%      *.1234229112,-.08579198697,.01321366966,.8970494003,9.136186247,   %! ASSUMED IN THE BIOT-SAVART INTEGRAL
%      *-38.19301215,21.73775846,-410.0783424,-69.90832690,-848.8543440,
%      *1.243288286,.2071721360,.05030555417,7.471332374,3.180533613,
%      *1.376743507,.1568504222,.02092910682,1.985148197,.3157139940,
%      *1.056309517,.1701395257,.1019870070,6.293740981,5.671824276,
%      *.1280772299,.02189060799,.01040696080,.1648265607,.04701592613,
%      *.01526400086,12.88384229,3.361775101,23.44173897/
A1 = 12.74640393;
A2 = -7.516393516;
A3 = -5.476233865;
A4 = 3.212704645;
A5 = -59.10926169;
A6 = 46.62198189;
A7 = -.01644280062;
A8 = .1234229112;
A9 = -.08579198697;
A10 = .01321366966;
A11 = .8970494003;
A12 = 9.136186247;
A13 = -38.19301215;
A14 = 21.73775846;
A15 = -410.0783424;
A16 = -69.90832690;
A17 = -848.8543440;
XK1 = 1.243288286;
AL1 = .2071721360;
DAL1 = .05030555417;
B1 = 7.471332374;
BE1 = 3.180533613;
XK2 = 1.376743507;
AL2 = .1568504222;
DAL2 = .02092910682;
BE2 = 1.985148197;
XK3 = .3157139940;
XK4 = 1.056309517;
AL3 = .1701395257;
DAL3 = .1019870070;
B3 = 6.293740981;
BE3 = 5.671824276;
AL4 = .1280772299;
DAL4 = .02189060799;
DG1 = .01040696080;
AL5 = .1648265607;
DAL5 = .04701592613;
DG2 = .01526400086;
C1 = 12.88384229;
C2 = 3.361775101;
C3 = 23.44173897;

SINT2=SINT^2;
COST2=COST^2;
SC=SINT*COST;
ALPHA=SINT2/R;
GAMMA=COST/R^2;

[F,FA,FS] = T01_FFS(ALPHA,AL1,DAL1);
D1=F^XK1/((R/B1)^BE1+1.D0);
D2=D1*COST2;

[F,FA,FS] = T01_FFS(ALPHA,AL2,DAL2);
D3=FA^XK2/R^BE2;
D4=D3*COST2;

[F,FA,FS] = T01_FFS(ALPHA,AL3,DAL3);
D5=FS^XK3*ALPHA^XK4/((R/B3)^BE3+1.D0);
D6=D5*COST2;

[F,FA,FS] = T01_FFS(GAMMA,0.D0,DG1);
FCC=(1.D0+((ALPHA-AL4)/DAL4)^2);
D7 =1.D0/FCC*FS;
D8 =D7/FCC;
D9 =D8/FCC;
D10=D9/FCC;

ARG=1.D0+((ALPHA-AL5)/DAL5)^2;
D11=1.D0/ARG/(1.D0+(GAMMA/DG2)^2);
D12=D11/ARG;
D13=D12/ARG;
D14=D13/ARG;

D15=1.D0/(R^4+C1^2);
D16=COST2/(R^4+C2^2);
D17=COST2^2/(R^4+C3^2);
% C
BT_PRC_Q=A1*D1+A2*D2+A3*D3+A4*D4+A5*D5+A6*D6+A7*D7+A8*D8+A9*D9+ ...
    A10*D10+A11*D11+A12*D12+A13*D13+A14*D14+A15*D15+A16*D16+A17*D17;
% C
% end of function BT_PRC_Q

% c
% c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

function [F,FA,FS] = T01_FFS(A,A0,DA);
%       SUBROUTINE FFS(A,A0,DA,F,FA,FS)
%       IMPLICIT  REAL * 8  (A - H, O - Z)
SQ1=sqrt((A+A0)^2+DA^2);
SQ2=sqrt((A-A0)^2+DA^2);
FA=2.D0/(SQ1+SQ2);
F=FA*A;
FS=0.5D0*(SQ1+SQ2)/(SQ1*SQ2)*(1.D0-F*F);
% end of function FFS

% C
% C||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
% C
% C
function [BX,BY,BZ] = T01_RC_SHIELD (A,PS,X_SC,X,Y,Z)
%         SUBROUTINE RC_SHIELD (A,PS,X_SC,X,Y,Z,BX,BY,BZ)
% C
% C   COMPUTES THE COMPONENTS OF THE SHIELDING FIELD FOR THE RING CURRENT
% C       (EITHER PARTIAL OR AXISYMMETRICAL)
% C   INPUT:   A - AN ARRAY CONTAINING THE HARMONIC COEFFICIENTS AND NONLINEAR PARAMETERS
% C            PS - GEODIPOLE TILT ANGLE IN RADIANS
% C            X_SC - SCALING FACTOR ( X_SC&gt 1 AND X_SC&lt 1 CORRESPOND TO LARGER/SMALLER
% C                  RING CURRENT, RESP.)
% C            X,Y,Z - POSITION IN RE (GSM COORDS)
% C   OUTPUT:  BX,BY,BZ - SHIELDING FIELD COMPONENTS (GSM)
% C
%	 IMPLICIT  REAL * 8  (A - H, O - Z)
%         DIMENSION A(86)
% C
FAC_SC=(X_SC+1.D0)^3;
% C
CPS=cos(PS);
SPS=sin(PS);

S3PS=2.D0*CPS;
% C
PST1=PS*A(85);
PST2=PS*A(86);

ST1=sin(PST1);
CT1=cos(PST1);
ST2=sin(PST2);
CT2=cos(PST2);

X1=X*CT1-Z*ST1;
Z1=X*ST1+Z*CT1;
X2=X*CT2-Z*ST2;
Z2=X*ST2+Z*CT2;
% C
L=0;
GX=0.D0;
GY=0.D0;
GZ=0.D0;
% C
for M=1:2,
    %         DO 1 M=1,2     %!    M=1 IS FOR THE 1ST SUM ("PERP." SYMMETRY)
    % C                           AND M=2 IS FOR THE SECOND SUM ("PARALL." SYMMETRY)
    for I=1:3,
        %             DO 2 I=1,3
        P=A(72+I);
        Q=A(78+I);
        CYPI=cos(Y/P);
        CYQI=cos(Y/Q);
        SYPI=sin(Y/P);
        SYQI=sin(Y/Q);
        % C
        for K=1:3,
            %                DO 3 K=1,3
            R=A(75+K);
            S=A(81+K);
            SZRK=sin(Z1/R);
            CZSK=cos(Z2/S);
            CZRK=cos(Z1/R);
            SZSK=sin(Z2/S);
            SQPR=sqrt(1.D0/P^2+1.D0/R^2);
            SQQS=sqrt(1.D0/Q^2+1.D0/S^2);
            EPR=exp(X1*SQPR);
            EQS=exp(X2*SQQS);
            % C
            for N=1:2,
                %                  DO 4 N=1,2  %! N=1 IS FOR THE FIRST PART OF EACH COEFFICIENT
                % C                                AND N=2 IS FOR THE SECOND ONE
                for NN=1:2,
                    %                    DO 5 NN=1,2 %!   NN = 1,2 FURTHER SPLITS THE COEFFICIENTS INTO 2 PARTS,
                    % C                                         TO TAKE INTO ACCOUNT THE SCALE FACTOR DEPENDENCE
                    
                    if (M==1) ,
                        FX=-SQPR*EPR*CYPI*SZRK  *FAC_SC;
                        FY=EPR*SYPI*SZRK/P   *FAC_SC;
                        FZ=-EPR*CYPI*CZRK/R  *FAC_SC;
                        if (N==1) ,
                            if (NN==1) ,
                                HX=FX;
                                HY=FY;
                                HZ=FZ;
                            else
                                HX=FX*X_SC;
                                HY=FY*X_SC;
                                HZ=FZ*X_SC;
                            end
                        else
                            if (NN==1) ,
                                HX=FX*CPS;
                                HY=FY*CPS;
                                HZ=FZ*CPS;
                            else
                                HX=FX*CPS*X_SC;
                                HY=FY*CPS*X_SC;
                                HZ=FZ*CPS*X_SC;
                            end
                        end
                        
                    else                            %!   M==2
                        FX=-SPS*SQQS*EQS*CYQI*CZSK  *FAC_SC;
                        FY=SPS/Q*EQS*SYQI*CZSK   *FAC_SC;
                        FZ=SPS/S*EQS*CYQI*SZSK   *FAC_SC;
                        if (N==1) ,
                            if (NN==1) ,
                                HX=FX;
                                HY=FY;
                                HZ=FZ;
                            else
                                HX=FX*X_SC;
                                HY=FY*X_SC;
                                HZ=FZ*X_SC;
                            end
                        else
                            if (NN==1) ,
                                HX=FX*S3PS;
                                HY=FY*S3PS;
                                HZ=FZ*S3PS;
                            else
                                HX=FX*S3PS*X_SC;
                                HY=FY*S3PS*X_SC;
                                HZ=FZ*S3PS*X_SC;
                            end
                        end
                    end
                    L=L+1;
                    
                    if (M==1) ,
                        HXR=HX*CT1+HZ*ST1;
                        HZR=-HX*ST1+HZ*CT1;
                    else
                        HXR=HX*CT2+HZ*ST2;
                        HZR=-HX*ST2+HZ*CT2;
                    end
                    
                    GX=GX+HXR*A(L);
                    GY=GY+HY *A(L);
                    GZ=GZ+HZR*A(L);
                end % 5
                
            end %  4   CONTINUE
        end %  3   CONTINUE
    end %  2   CONTINUE
end %  1   CONTINUE

BX=GX;
BY=GY;
BZ=GZ;

% end of function RC_SHIELD

% C
% c===========================================================================
% c
function [BX,BY,BZ] = T01_DIPOLE (PS,X,Y,Z)
%      SUBROUTINE DIPOLE (PS,X,Y,Z,BX,BY,BZ)
% C
% C     THIS IS A DOUBLE PRECISION ROUTINE, OTHERWISE IDENTICAL TO THE S/R DIP OF GEOPACK
% C
% C     CALCULATES GSM COMPONENTS OF A GEODIPOLE FIELD WITH THE DIPOLE MOMENT
% C     CORRESPONDING TO THE EPOCH OF 2000.
% C
% C------INPUT PARAMETERS:
% C       PS - GEODIPOLE TILT ANGLE IN RADIANS,
% C       X,Y,Z - GSM COORDINATES IN RE (1 RE = 6371.2 km)
% C
% C----OUTPUT PARAMETERS:
% C     BX,BY,BZ - FIELD COMPONENTS IN GSM SYSTEM, IN NANOTESLA.
% C
% C    LAST MODIFICATION: JAN. 5, 2001. THE VALUE OF THE DIPOLE MOMENT WAS UPDATED TO 2000.
% C      AND A "SAVE" STATEMENT HAS BEEN ADDED, TO AVOID POTENTIAL PROBLEMS WITH SOME
% C      FORTRAN COMPILERS
% C
% C    WRITTEN BY: N. A. TSYGANENKO
% C
%      IMPLICIT REAL*8 (A-H,O-Z)
persistent M PSI SPS CPS
if isempty(M),
    M = 0;
end
if isempty(PSI),
    PSI = 5.D0;
end
if~((M==1) & (abs(PS-PSI)<1.D-5)) ,    %!   THIS IS TO AVOID MULTIPLE CALCULATIONS
    SPS=sin(PS);                                  %!   OF SIN(PS) AND COS(PS), IF THE ANGLE PS
    CPS=cos(PS);                                  %!   REMAINS UNCHANGED
    PSI=PS;
    M=1;
end
P=X^2;
U=Z^2;
V=3.D0*Z*X;
T=Y^2;
Q=30115.D0/sqrt(P+T+U)^5;
BX=Q*((T+U-2.D0*P)*SPS-V*CPS);
BY=-3.D0*Y*Q*(X*SPS+Z*CPS);
BZ=Q*((P+T-2.D0*U)*CPS-V*SPS);

% end of function DIPOLE


% c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@