117c117
<       CALL EXTALL (0,0,0,0,A,43,PDYN,DST_AST,BYIMF,BZIMF,G1,G2,
---
>       CALL T01EXTALL (0,0,0,0,A,43,PDYN,DST_AST,BYIMF,BZIMF,G1,G2,
131c131
<       SUBROUTINE EXTALL (IOPGEN,IOPT,IOPB,IOPR,A,NTOT,
---
>       SUBROUTINE T01EXTALL (IOPGEN,IOPT,IOPB,IOPR,A,NTOT,
165c165
<       COMMON /RH0/ RH0
---
>       COMMON /RH0T01/ RH0T01
168c168
<       DATA DSIG /0.003D0/, RH0,RH2 /8.0D0,-5.2D0/
---
>       DATA DSIG /0.003D0/, RH0T01,RH2 /8.0D0,-5.2D0/
171c171
<       RH0=A(40)
---
>       RH0T01=A(40)
218c218
<   1   XSOLD=XSS      !   BEGIN ITERATIVE SEARCH OF UNWARPED COORDS (TO FIND SIGMA)
---
>   1   XSOLD=XSS      !   BEGIN ITERATIVE SEARCH OF UNT01WARPED COORDS (TO FIND SIGMA)
221c221
<       RH=RH0+RH2*(ZSS/R)**2
---
>       RH=RH0T01+RH2*(ZSS/R)**2
249c249
<          CALL SHLCAR3X301(XX,YY,ZZ,PS,CFX,CFY,CFZ)         !  DIPOLE SHIELDING FIELD
---
>          CALL T01SHLCAR3X301(XX,YY,ZZ,PS,CFX,CFY,CFZ)         !  DIPOLE SHIELDING FIELD
264c264
<          CALL DEFORMED (IOPT,PS,XX,YY,ZZ,                !  TAIL FIELD (THREE MODES)
---
>          CALL T01DEFORMED (IOPT,PS,XX,YY,ZZ,                !  TAIL FIELD (THREE MODES)
278c278
<          CALL BIRK_TOT (IOPB,PS,XX,YY,ZZ,BXR11,BYR11,BZR11,BXR12,BYR12,
---
>          CALL T01BIRK_TOT (IOPB,PS,XX,YY,ZZ,BXR11,BYR11,BZR11,BXR12,BYR12,
301c301
<          CALL FULL_RC(IOPR,PS,XX,YY,ZZ,BXSRC,BYSRC,BZSRC,BXPRC,BYPRC,
---
>          CALL T01FULL_RC(IOPR,PS,XX,YY,ZZ,BXSRC,BYSRC,BZSRC,BXPRC,BYPRC,
369c369
<        CALL DIPOLE01 (PS,X,Y,Z,QX,QY,QZ)
---
>        CALL T01DIPOLE01 (PS,X,Y,Z,QX,QY,QZ)
378c378
<                 CALL DIPOLE01 (PS,X,Y,Z,QX,QY,QZ)
---
>                 CALL T01DIPOLE01 (PS,X,Y,Z,QX,QY,QZ)
388c388
<          SUBROUTINE  SHLCAR3X301(X,Y,Z,PS,BX,BY,BZ)
---
>          SUBROUTINE  T01SHLCAR3X301(X,Y,Z,PS,BX,BY,BZ)
714c714
< c
---
> 
718c718
<       SUBROUTINE DEFORMED (IOPT,PS,X,Y,Z,BX1,BY1,BZ1,BX2,BY2,BZ2)
---
>       SUBROUTINE T01DEFORMED (IOPT,PS,X,Y,Z,BX1,BY1,BZ1,BX2,BY2,BZ2)
726c726
< C    WARPING IN Y-Z (DONE BY THE S/R WARPED) AND BENDING IN X-Z (DONE BY THIS SUBROUTINE)
---
> C    WARPING IN Y-Z (DONE BY THE S/R T01WARPED) AND BENDING IN X-Z (DONE BY THIS SUBROUTINE)
729c729
<       COMMON /RH0/ RH0
---
>       COMMON /RH0T01/ RH0T01
732c732
< C  RH0,RH1,RH2, AND IEPS CONTROL THE TILT-RELATED DEFORMATION OF THE TAIL FIELD
---
> C  RH0T01,RH1,RH2, AND IEPS CONTROL THE TILT-RELATED DEFORMATION OF THE TAIL FIELD
739c739
<       RH=RH0+RH2*ZR**2
---
>       RH=RH0T01+RH2*ZR**2
771c771,772
<       CALL WARPED(IOPT,PS,XAS,Y,ZAS,BXAS1,BYAS1,BZAS1,BXAS2,BYAS2,BZAS2)
---
>       CALL T01WARPED(IOPT,PS,XAS,Y,ZAS,BXAS1,BYAS1,BZAS1,BXAS2,BYAS2,
>      + BZAS2)
783c784
< C
---
> cC
786c787
<       SUBROUTINE WARPED (IOPT,PS,X,Y,Z,BX1,BY1,BZ1,BX2,BY2,BZ2)
---
>       SUBROUTINE T01WARPED (IOPT,PS,X,Y,Z,BX1,BY1,BZ1,BX2,BY2,BZ2)
788,790c789,791
< 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   CALCULATES GSM COMPONENTS OF THE T01WARPED FIELD FOR TWO TAIL UNIT MODES.
> C   THE WARPING DEFORMATION IS IMPOSED ON THE UNT01WARPED FIELD, COMPUTED
> C   BY THE S/R "UNT01WARPED".  THE WARPING PARAMETER G WAS OBTAINED BY LEAST
830c831
<       CALL UNWARPED (IOPT,X,YAS,ZAS,BX_AS1,BY_AS1,BZ_AS1,
---
>       CALL UNT01WARPED (IOPT,X,YAS,ZAS,BX_AS1,BY_AS1,BZ_AS1,
858c859
<        SUBROUTINE UNWARPED (IOPT,X,Y,Z,BX1,BY1,BZ1,BX2,BY2,BZ2)
---
>        SUBROUTINE UNT01WARPED (IOPT,X,Y,Z,BX1,BY1,BZ1,BX2,BY2,BZ2)
914,915c915,917
<       CALL TAILDISK01(D0SC1,DELTADX1,DELTADY,XSC1,YSC1,ZSC1,FX1,FY1,FZ1)
<       CALL SHLCAR5X5(A1,X,Y,Z,DXSHIFT1,HX1,HY1,HZ1)
---
>       CALL T01TAILDISK01(D0SC1,DELTADX1,DELTADY,XSC1,YSC1,ZSC1,FX1,FY1,
>      + FZ1)
>       CALL T01SHLCAR5X5(A1,X,Y,Z,DXSHIFT1,HX1,HY1,HZ1)
933,934c935,937
<       CALL TAILDISK01(D0SC2,DELTADX2,DELTADY,XSC2,YSC2,ZSC2,FX2,FY2,FZ2)
<       CALL SHLCAR5X5(A2,X,Y,Z,DXSHIFT2,HX2,HY2,HZ2)
---
>       CALL T01TAILDISK01(D0SC2,DELTADX2,DELTADY,XSC2,YSC2,ZSC2,FX2,FY2,
>      + FZ2)
>       CALL T01SHLCAR5X5(A2,X,Y,Z,DXSHIFT2,HX2,HY2,HZ2)
952c955
<       SUBROUTINE TAILDISK01(D0,DELTADX,DELTADY,X,Y,Z,BX,BY,BZ)
---
>       SUBROUTINE T01TAILDISK01(D0,DELTADX,DELTADY,X,Y,Z,BX,BY,BZ)
1043c1046
<          SUBROUTINE  SHLCAR5X5(A,X,Y,Z,DSHIFT,HX,HY,HZ)
---
>          SUBROUTINE  T01SHLCAR5X5(A,X,Y,Z,DSHIFT,HX,HY,HZ)
1097c1100
<       SUBROUTINE BIRK_TOT (IOPB,PS,X,Y,Z,BX11,BY11,BZ11,BX12,BY12,BZ12,
---
>       SUBROUTINE T01BIRK_TOT (IOPB,PS,X,Y,Z,BX11,BY11,BZ11,BX12,BY12,BZ12,
1107c1110
<       COMMON /BIRKPAR/ XKAPPA1,XKAPPA2   !  INPUT PARAMETERS, SPECIFIED FROM S/R EXTALL
---
>       COMMON /BIRKPAR/ XKAPPA1,XKAPPA2   !  INPUT PARAMETERS, SPECIFIED FROM S/R T01EXTALL
1186,1187c1189,1190
<       XKAPPA=XKAPPA1        !  FORWARDED IN BIRK_1N2
<       X_SC=XKAPPA1-1.1D0    !  FORWARDED IN BIRK_SHL
---
>       XKAPPA=XKAPPA1        !  FORWARDED IN T01BIRK_1N2
>       X_SC=XKAPPA1-1.1D0    !  FORWARDED IN T01BIRK_SHL
1191,1192c1194,1195
<       CALL BIRK_1N2 (1,1,PS,X,Y,Z,FX11,FY11,FZ11)           !  REGION 1, MODE 1
<       CALL BIRK_SHL (SH11,PS,X_SC,X,Y,Z,HX11,HY11,HZ11)
---
>       CALL T01BIRK_1N2 (1,1,PS,X,Y,Z,FX11,FY11,FZ11)           !  REGION 1, MODE 1
>       CALL T01BIRK_SHL (SH11,PS,X_SC,X,Y,Z,HX11,HY11,HZ11)
1197,1198c1200,1201
<       CALL BIRK_1N2 (1,2,PS,X,Y,Z,FX12,FY12,FZ12)           !  REGION 1, MODE 2
<       CALL BIRK_SHL (SH12,PS,X_SC,X,Y,Z,HX12,HY12,HZ12)
---
>       CALL T01BIRK_1N2 (1,2,PS,X,Y,Z,FX12,FY12,FZ12)           !  REGION 1, MODE 2
>       CALL T01BIRK_SHL (SH12,PS,X_SC,X,Y,Z,HX12,HY12,HZ12)
1205,1206c1208,1209
<       XKAPPA=XKAPPA2        !  FORWARDED IN BIRK_1N2
<       X_SC=XKAPPA2-1.0D0    !  FORWARDED IN BIRK_SHL
---
>       XKAPPA=XKAPPA2        !  FORWARDED IN T01BIRK_1N2
>       X_SC=XKAPPA2-1.0D0    !  FORWARDED IN T01BIRK_SHL
1210,1211c1213,1214
<       CALL BIRK_1N2 (2,1,PS,X,Y,Z,FX21,FY21,FZ21)           !  REGION 2, MODE 1
<       CALL BIRK_SHL (SH21,PS,X_SC,X,Y,Z,HX21,HY21,HZ21)
---
>       CALL T01BIRK_1N2 (2,1,PS,X,Y,Z,FX21,FY21,FZ21)           !  REGION 2, MODE 1
>       CALL T01BIRK_SHL (SH21,PS,X_SC,X,Y,Z,HX21,HY21,HZ21)
1216,1217c1219,1220
<       CALL BIRK_1N2 (2,2,PS,X,Y,Z,FX22,FY22,FZ22)           !  REGION 2, MODE 2
<       CALL BIRK_SHL (SH22,PS,X_SC,X,Y,Z,HX22,HY22,HZ22)
---
>       CALL T01BIRK_1N2 (2,2,PS,X,Y,Z,FX22,FY22,FZ22)           !  REGION 2, MODE 2
>       CALL T01BIRK_SHL (SH22,PS,X_SC,X,Y,Z,HX22,HY22,HZ22)
1230c1233
<       SUBROUTINE BIRK_1N2 (NUMB,MODE,PS,X,Y,Z,BX,BY,BZ)        !   NB# 6, P.60
---
>       SUBROUTINE T01BIRK_1N2 (NUMB,MODE,PS,X,Y,Z,BX,BY,BZ)        !   NB# 6, P.60
1341,1342c1344,1345
<         IF (MODE.EQ.1) CALL TWOCONES (A11,XS,Ysc,ZS,BXS,BYAS,BZS)
<         IF (MODE.EQ.2) CALL TWOCONES (A12,XS,Ysc,ZS,BXS,BYAS,BZS)
---
>         IF (MODE.EQ.1) CALL T01TWOCONES (A11,XS,Ysc,ZS,BXS,BYAS,BZS)
>         IF (MODE.EQ.2) CALL T01TWOCONES (A12,XS,Ysc,ZS,BXS,BYAS,BZS)
1344,1345c1347,1348
<         IF (MODE.EQ.1) CALL TWOCONES (A21,XS,Ysc,ZS,BXS,BYAS,BZS)
<         IF (MODE.EQ.2) CALL TWOCONES (A22,XS,Ysc,ZS,BXS,BYAS,BZS)
---
>         IF (MODE.EQ.1) CALL T01TWOCONES (A21,XS,Ysc,ZS,BXS,BYAS,BZS)
>         IF (MODE.EQ.2) CALL T01TWOCONES (A22,XS,Ysc,ZS,BXS,BYAS,BZS)
1364c1367
<       SUBROUTINE TWOCONES (A,X,Y,Z,BX,BY,BZ)
---
>       SUBROUTINE T01TWOCONES (A,X,Y,Z,BX,BY,BZ)
1372,1373c1375,1376
<       CALL ONE_CONE (A,X,Y,Z,BXN,BYN,BZN)
<       CALL ONE_CONE (A,X,-Y,-Z,BXS,BYS,BZS)
---
>       CALL T01ONE_CONE (A,X,Y,Z,BXN,BYN,BZN)
>       CALL T01ONE_CONE (A,X,-Y,-Z,BXS,BYS,BZS)
1383c1386
<       SUBROUTINE ONE_CONE(A,X,Y,Z,BX,BY,BZ)
---
>       SUBROUTINE T01ONE_CONE(A,X,Y,Z,BX,BY,BZ)
1385c1388
< c  RETURNS FIELD COMPONENTS FOR A DEFORMED CONICAL CURRENT SYSTEM, FITTED TO A BIOSAVART FIELD
---
> c  RETURNS FIELD COMPONENTS FOR A T01DEFORMED CONICAL CURRENT SYSTEM, FITTED TO A BIOSAVART FIELD
1412c1415
<        CALL FIALCOS (RS,THETAS,PHIS,BTAST,BFAST,M,THETA0,DTHETA)    !   MODE #M
---
>        CALL T01FIALCOS (RS,THETAS,PHIS,BTAST,BFAST,M,THETA0,DTHETA)    !   MODE #M
1444,1472c1447,1475
< C=====================================================================================
<       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/DSQRT(R**2+A(11)**2)+A(4)*R/(R**2+A(12)**2)
<      *+(A(5)+A(6)/R+A(7)*R/DSQRT(R**2+A(13)**2)+A(8)*R/(R**2+A(14)**2))*
<      * DCOS(THETA)
<      *+(A(9)*R/DSQRT(R**2+A(15)**2)+A(10)*R/(R**2+A(16)**2)**2)
<      * *DCOS(2.D0*THETA)
< C
<       RETURN
<       END
< C
< C-----------------------------------------------------------------------------
< C
<       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/DSQRT(R**2+A(27)**2))*DSIN(THETA)
<      * +(A(21)+A(22)*R/DSQRT(R**2+A(28)**2)
<      *                +A(23)*R/(R**2+A(29)**2))*DSIN(2.D0*THETA)
<      * +(A(24)+A(25)/R+A(26)*R/(R**2+A(30)**2))*DSIN(3.D0*THETA)
< C
<       RETURN
<       END
< C
---
> cC=====================================================================================
> c      DOUBLE PRECISION FUNCTION R_S(A,R,THETA)
> c      IMPLICIT REAL*8 (A-H,O-Z)
> c      DIMENSION A(31)
> cC
> c      R_S=R+A(2)/R+A(3)*R/DSQRT(R**2+A(11)**2)+A(4)*R/(R**2+A(12)**2)
> c     *+(A(5)+A(6)/R+A(7)*R/DSQRT(R**2+A(13)**2)+A(8)*R/(R**2+A(14)**2))*
> c     * DCOS(THETA)
> c     *+(A(9)*R/DSQRT(R**2+A(15)**2)+A(10)*R/(R**2+A(16)**2)**2)
> c     * *DCOS(2.D0*THETA)
> cC
> c      RETURN
> c      END
> cC
> cC-----------------------------------------------------------------------------
> cC
> c      DOUBLE PRECISION FUNCTION THETA_S(A,R,THETA)
> c      IMPLICIT REAL*8 (A-H,O-Z)
> c      DIMENSION A(31)
> cc
> c      THETA_S=THETA+(A(17)+A(18)/R+A(19)/R**2
> c     *                +A(20)*R/DSQRT(R**2+A(27)**2))*DSIN(THETA)
> c     * +(A(21)+A(22)*R/DSQRT(R**2+A(28)**2)
> c     *                +A(23)*R/(R**2+A(29)**2))*DSIN(2.D0*THETA)
> c     * +(A(24)+A(25)/R+A(26)*R/(R**2+A(30)**2))*DSIN(3.D0*THETA)
> cC
> c      RETURN
> c      END
> cC
1475c1478
<       SUBROUTINE FIALCOS(R,THETA,PHI,BTHETA,BPHI,N,THETA0,DT)
---
>       SUBROUTINE T01FIALCOS(R,THETA,PHI,BTHETA,BPHI,N,THETA0,DT)
1554c1557
<          SUBROUTINE BIRK_SHL (A,PS,X_SC,X,Y,Z,BX,BY,BZ)
---
>          SUBROUTINE T01BIRK_SHL (A,PS,X_SC,X,Y,Z,BX,BY,BZ)
1691c1694
<       SUBROUTINE FULL_RC (IOPR,PS,X,Y,Z,BXSRC,BYSRC,BZSRC,BXPRC,BYPRC,
---
>       SUBROUTINE T01FULL_RC (IOPR,PS,X,Y,Z,BXSRC,BYSRC,BZSRC,BXPRC,BYPRC,
1751c1754
<         CALL SRC_PRC (IOPR,SC_SY,SC_PR,PHI,PS,X,Y,Z,HXSRC,HYSRC,HZSRC,
---
>         CALL T01SRC_PRC (IOPR,SC_SY,SC_PR,PHI,PS,X,Y,Z,HXSRC,HYSRC,HZSRC,
1756c1759
<           CALL RC_SHIELD (C_SY,PS,X_SC,X,Y,Z,FSX,FSY,FSZ)
---
>           CALL T01RC_SHIELD (C_SY,PS,X_SC,X,Y,Z,FSX,FSY,FSZ)
1765c1768
<           CALL RC_SHIELD (C_PR,PS,X_SC,X,Y,Z,FPX,FPY,FPZ)
---
>           CALL T01RC_SHIELD (C_PR,PS,X_SC,X,Y,Z,FPX,FPY,FPZ)
1784c1787
<        SUBROUTINE SRC_PRC (IOPR,SC_SY,SC_PR,PHI,PS,X,Y,Z,BXSRC,BYSRC,
---
>        SUBROUTINE T01SRC_PRC (IOPR,SC_SY,SC_PR,PHI,PS,X,Y,Z,BXSRC,BYSRC,
1838c1841
<         IF (IOPR.LE.1) CALL RC_SYMM(XTS,YTS,ZTS,BXS,BYS,BZS)
---
>         IF (IOPR.LE.1) CALL T01RC_SYMM(XTS,YTS,ZTS,BXS,BYS,BZS)
1840c1843
<      *                 CALL PRC_SYMM(XTA,YTA,ZTA,BXA_S,BYA_S,BZA_S)
---
>      *                 CALL PT01RC_SYMM(XTA,YTA,ZTA,BXA_S,BYA_S,BZA_S)
1851c1854
<      *                 CALL PRC_QUAD(XR,YR,ZTA,BXA_QR,BYA_QR,BZA_Q)
---
>      *                 CALL T01PRC_QUAD(XR,YR,ZTA,BXA_QR,BYA_QR,BZA_Q)
1879c1882
<       SUBROUTINE RC_SYMM (X,Y,Z,BX,BY,BZ)
---
>       SUBROUTINE T01RC_SYMM (X,Y,Z,BX,BY,BZ)
1921,2036c1924,2039
< c
< c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
< C
<       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/
< 
<       PROX=.FALSE.
<       SINT1=SINT
<       COST1=COST
<       IF (SINT1.LT.1.D-2) THEN  !  TOO CLOSE TO Z-AXIS;  USE LINEAR INTERPOLATION BETWEEN SINT=0 & SINT=0.01
<         SINT1=1.D-2
<         COST1=.99994999875
<         PROX=.TRUE.
<       ENDIF
< 
<          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.LT.-500.D0) THEN        !   TO PREVENT "FLOATING UNDERFLOW" CRASHES
<            DEXP1=0.D0
<          ELSE
<            DEXP1=DEXP(ARG1)
<          ENDIF
< 
<          IF (ARG2.LT.-500.D0) THEN
<            DEXP2=0.D0
<          ELSE
<            DEXP2=DEXP(ARG2)
<          ENDIF
< 
<          IF (ARG3.LT.-500.D0) THEN
<            DEXP3=0.D0
<          ELSE
<            DEXP3=DEXP(ARG3)
<          ENDIF
< 
< 
<          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=(DSQRT(F)+ALSQH)**(1.D0/3.D0)
<          C=Q-4.D0*GAMMAS2**(1.D0/3.D0)/(3.D0*Q)
<          IF (C.LT.0.D0) C=0.D0
<          G=DSQRT(C**2+4.D0*GAMMAS2**(1.D0/3.D0))
<          RS=4.D0/((DSQRT(2.D0*G-C)+DSQRT(C))*(G+C))
<          COSTS=GAMMA_S*RS**2
<          SINTS=DSQRT(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=DLOG(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=DLOG(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   !   LINEAR INTERPOLATION, IF TOO CLOSE TO THE Z-AXIS
< C
<        RETURN
<        END
---
> cc
> cc&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
> cC
> c      DOUBLE PRECISION FUNCTION AP(R,SINT,COST)
> cC
> cC      Calculates azimuthal component of the vector potential of the symmetric
> cc  part of the model ring current.
> cC
> c      IMPLICIT  REAL * 8  (A - H, O - Z)
> c      LOGICAL PROX   !  INDICATES WHETHER WE ARE TOO CLOSE TO THE AXIS OF SYMMETRY, WHERE THE INVERSION
> cC                                                             OF DIPOLAR COORDINATES BECOMES INACCURATE
> c      DATA A1,A2,RRC1,DD1,RRC2,DD2,P1,R1,DR1,DLA1,P2,R2,DR2,DLA2,P3,
> c     *R3,DR3/-563.3722359,425.0891691,4.150588549,2.266150226,
> c     * 3.334503403,3.079071195,.02602428295,8.937790598,3.327934895,
> c     *.4487061833,.09125832351,6.243029867,1.750145910,.4181957162,
> c     *.06106691992,2.079908581,.6828548533/
> 
> c      PROX=.FALSE.
> c      SINT1=SINT
> c      COST1=COST
> c      IF (SINT1.LT.1.D-2) THEN  !  TOO CLOSE TO Z-AXIS;  USE LINEAR INTERPOLATION BETWEEN SINT=0 & SINT=0.01
> c        SINT1=1.D-2
> c        COST1=.99994999875
> c        PROX=.TRUE.
> c      ENDIF
> 
> c         ALPHA=SINT1**2/R         !  R,THETA -> ALPHA,GAMMA
> c         GAMMA=COST1/R**2
> 
> c         ARG1=-((R-R1)/DR1)**2-(COST1/DLA1)**2
> c         ARG2=-((R-R2)/DR2)**2-(COST1/DLA2)**2
> c         ARG3=-((R-R3)/DR3)**2
> 
> c         IF (ARG1.LT.-500.D0) THEN        !   TO PREVENT "FLOATING UNDERFLOW" CRASHES
> c           DEXP1=0.D0
> c         ELSE
> c           DEXP1=DEXP(ARG1)
> c         ENDIF
> 
> c         IF (ARG2.LT.-500.D0) THEN
> c           DEXP2=0.D0
> c         ELSE
> c           DEXP2=DEXP(ARG2)
> c         ENDIF
> 
> c         IF (ARG3.LT.-500.D0) THEN
> c           DEXP3=0.D0
> c         ELSE
> c           DEXP3=DEXP(ARG3)
> c         ENDIF
> 
> 
> c         ALPHA_S=ALPHA*(1.D0+P1*DEXP1+P2*DEXP2+P3*DEXP3)     !  ALPHA -> ALPHA_S  (T01DEFORMED)
> 
> c         GAMMA_S=GAMMA
> c         GAMMAS2=GAMMA_S**2
> 
> 
> c         ALSQH=ALPHA_S**2/2.D0            !  ALPHA_S,GAMMA_S -> RS,SINTS,COSTS
> c         F=64.D0/27.D0*GAMMAS2+ALSQH**2
> c         Q=(DSQRT(F)+ALSQH)**(1.D0/3.D0)
> c         C=Q-4.D0*GAMMAS2**(1.D0/3.D0)/(3.D0*Q)
> c         IF (C.LT.0.D0) C=0.D0
> c         G=DSQRT(C**2+4.D0*GAMMAS2**(1.D0/3.D0))
> c         RS=4.D0/((DSQRT(2.D0*G-C)+DSQRT(C))*(G+C))
> c         COSTS=GAMMA_S*RS**2
> c         SINTS=DSQRT(1.D0-COSTS**2)
> c         RHOS=RS*SINTS
> c         RHOS2=RHOS**2
> c         ZS=RS*COSTS
> cC
> cc  1st loop:
> 
> c         P=(RRC1+RHOS)**2+ZS**2+DD1**2
> c         XK2=4.D0*RRC1*RHOS/P
> c         XK=SQRT(XK2)
> c         XKRHO12=XK*SQRT(RHOS)     !   SEE NB#4, P.3
> cC
> c      XK2S=1.D0-XK2
> c      DL=DLOG(1.D0/XK2S)
> c      ELK=1.38629436112d0+XK2S*(0.09666344259D0+XK2S*(0.03590092383+
> c     *     XK2S*(0.03742563713+XK2S*0.01451196212))) +DL*
> c     *     (0.5D0+XK2S*(0.12498593597D0+XK2S*(0.06880248576D0+
> c     *     XK2S*(0.03328355346D0+XK2S*0.00441787012D0))))
> c      ELE=1.D0+XK2S*(0.44325141463D0+XK2S*(0.0626060122D0+XK2S*
> c     *      (0.04757383546D0+XK2S*0.01736506451D0))) +DL*
> c     *     XK2S*(0.2499836831D0+XK2S*(0.09200180037D0+XK2S*
> c     *       (0.04069697526D0+XK2S*0.00526449639D0)))
> cC
> c      APHI1=((1.D0-XK2*0.5D0)*ELK-ELE)/XKRHO12
> cc
> cc  2nd loop:
> 
> c         P=(RRC2+RHOS)**2+ZS**2+DD2**2
> c         XK2=4.D0*RRC2*RHOS/P
> c         XK=SQRT(XK2)
> c         XKRHO12=XK*SQRT(RHOS)     !   SEE NB#4, P.3
> cC
> c      XK2S=1.D0-XK2
> c      DL=DLOG(1.D0/XK2S)
> c      ELK=1.38629436112d0+XK2S*(0.09666344259D0+XK2S*(0.03590092383+
> c     *     XK2S*(0.03742563713+XK2S*0.01451196212))) +DL*
> c     *     (0.5D0+XK2S*(0.12498593597D0+XK2S*(0.06880248576D0+
> c     *     XK2S*(0.03328355346D0+XK2S*0.00441787012D0))))
> c      ELE=1.D0+XK2S*(0.44325141463D0+XK2S*(0.0626060122D0+XK2S*
> c     *      (0.04757383546D0+XK2S*0.01736506451D0))) +DL*
> c     *     XK2S*(0.2499836831D0+XK2S*(0.09200180037D0+XK2S*
> c     *       (0.04069697526D0+XK2S*0.00526449639D0)))
> cC
> c       APHI2=((1.D0-XK2*0.5D0)*ELK-ELE)/XKRHO12
> 
> c       AP=A1*APHI1+A2*APHI2
> c       IF (PROX) AP=AP*SINT/SINT1   !   LINEAR INTERPOLATION, IF TOO CLOSE TO THE Z-AXIS
> cC
> c       RETURN
> c       END
2040c2043
<       SUBROUTINE PRC_SYMM (X,Y,Z,BX,BY,BZ)
---
>       SUBROUTINE PT01RC_SYMM (X,Y,Z,BX,BY,BZ)
2085,2188c2088,2191
< C
<       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/
< 
<       PROX=.FALSE.
<       SINT1=SINT
<       COST1=COST
<       IF (SINT1.LT.1.D-2) THEN  !  TOO CLOSE TO Z-AXIS;  USE LINEAR INTERPOLATION BETWEEN SINT=0 & SINT=0.01
<         SINT1=1.D-2
<         COST1=.99994999875
<         PROX=.TRUE.
<       ENDIF
< 
<          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
<      * *DEXP(-(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)
<      * *DEXP(-((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=(DSQRT(F)+ALSQH)**(1.D0/3.D0)
<          C=Q-4.D0*GAMMAS2**(1.D0/3.D0)/(3.D0*Q)
<          IF (C.LT.0.D0) C=0.D0
<          G=DSQRT(C**2+4.D0*GAMMAS2**(1.D0/3.D0))
<          RS=4.D0/((DSQRT(2.D0*G-C)+DSQRT(C))*(G+C))
<          COSTS=GAMMA_S*RS**2
<          SINTS=DSQRT(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=DLOG(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=DLOG(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   !   LINEAR INTERPOLATION, IF TOO CLOSE TO THE Z-AXIS
< C
<       RETURN
<       END
< C
---
> cC
> c      DOUBLE PRECISION FUNCTION APPRC(R,SINT,COST)
> cC
> cC      Calculates azimuthal component of the vector potential of the symmetric
> cc  part of the model PARTIAL ring current.
> cC
> c      IMPLICIT  REAL * 8  (A - H, O - Z)
> c      LOGICAL PROX
> c      DATA A1,A2,RRC1,DD1,RRC2,DD2,P1,ALPHA1,DAL1,BETA1,DG1,P2,ALPHA2,
> c     * DAL2,BETA2,DG2,BETA3,P3,ALPHA3,DAL3,BETA4,DG3,BETA5,Q0,Q1,ALPHA4,
> c     * DAL4,DG4,Q2,ALPHA5,DAL5,DG5,BETA6,BETA7
> c     * /-80.11202281,12.58246758,6.560486035,1.930711037,3.827208119,
> c     *.7789990504,.3058309043,.1817139853,.1257532909,3.422509402,
> c     *.04742939676,-4.800458958,-.02845643596,.2188114228,2.545944574,
> c     *.00813272793,.35868244,103.1601001,-.00764731187,.1046487459,
> c     *2.958863546,.01172314188,.4382872938,.01134908150,14.51339943,
> c     *.2647095287,.07091230197,.01512963586,6.861329631,.1677400816,
> c     *.04433648846,.05553741389,.7665599464,.7277854652/
> 
> c      PROX=.FALSE.
> c      SINT1=SINT
> c      COST1=COST
> c      IF (SINT1.LT.1.D-2) THEN  !  TOO CLOSE TO Z-AXIS;  USE LINEAR INTERPOLATION BETWEEN SINT=0 & SINT=0.01
> c        SINT1=1.D-2
> c        COST1=.99994999875
> c        PROX=.TRUE.
> c      ENDIF
> 
> c         ALPHA=SINT1**2/R         !  R,THETA -> ALPHA,GAMMA
> c         GAMMA=COST1/R**2
> 
> c         ALPHA_S=ALPHA*(1.D0+P1/(1.D0+((ALPHA-ALPHA1)/DAL1)**2)**BETA1
> c     * *DEXP(-(GAMMA/DG1)**2)
> c     *+P2*(ALPHA-ALPHA2)/(1.D0+((ALPHA-ALPHA2)/DAL2)**2)**BETA2
> c     */(1.D0+(GAMMA/DG2)**2)**BETA3
> c     *+P3*(ALPHA-ALPHA3)**2/(1.D0+((ALPHA-ALPHA3)/DAL3)**2)**BETA4
> c     */(1.D0+(GAMMA/DG3)**2)**BETA5)     !  ALPHA -> ALPHA_S  (T01DEFORMED)
> 
> c         GAMMA_S=GAMMA*(1.D0+Q0+Q1*(ALPHA-ALPHA4)
> c     * *DEXP(-((ALPHA-ALPHA4)/DAL4)**2-(GAMMA/DG4)**2)   !  GAMMA -> GAMMA_  (T01DEFORMED)
> c     * +Q2*(ALPHA-ALPHA5)/(1.D0+((ALPHA-ALPHA5)/DAL5)**2)**BETA6
> c     * /(1.D0+(GAMMA/DG5)**2)**BETA7)
> 
> c         GAMMAS2=GAMMA_S**2
> 
> c         ALSQH=ALPHA_S**2/2.D0                            !  ALPHA_S,GAMMA_S -> RS,SINTS,COSTS
> c         F=64.D0/27.D0*GAMMAS2+ALSQH**2
> c         Q=(DSQRT(F)+ALSQH)**(1.D0/3.D0)
> c         C=Q-4.D0*GAMMAS2**(1.D0/3.D0)/(3.D0*Q)
> c         IF (C.LT.0.D0) C=0.D0
> c         G=DSQRT(C**2+4.D0*GAMMAS2**(1.D0/3.D0))
> c         RS=4.D0/((DSQRT(2.D0*G-C)+DSQRT(C))*(G+C))
> c         COSTS=GAMMA_S*RS**2
> c         SINTS=DSQRT(1.D0-COSTS**2)
> c         RHOS=RS*SINTS
> c         RHOS2=RHOS**2
> c         ZS=RS*COSTS
> cC
> cc  1st loop:
> 
> c         P=(RRC1+RHOS)**2+ZS**2+DD1**2
> c         XK2=4.D0*RRC1*RHOS/P
> c         XK=SQRT(XK2)
> c         XKRHO12=XK*SQRT(RHOS)     !    NB#4, P.3
> cC
> c      XK2S=1.D0-XK2
> c      DL=DLOG(1.D0/XK2S)
> c      ELK=1.38629436112d0+XK2S*(0.09666344259D0+XK2S*(0.03590092383+
> c     *     XK2S*(0.03742563713+XK2S*0.01451196212))) +DL*
> c     *     (0.5D0+XK2S*(0.12498593597D0+XK2S*(0.06880248576D0+
> c     *     XK2S*(0.03328355346D0+XK2S*0.00441787012D0))))
> c      ELE=1.D0+XK2S*(0.44325141463D0+XK2S*(0.0626060122D0+XK2S*
> c     *      (0.04757383546D0+XK2S*0.01736506451D0))) +DL*
> c     *     XK2S*(0.2499836831D0+XK2S*(0.09200180037D0+XK2S*
> c     *       (0.04069697526D0+XK2S*0.00526449639D0)))
> cC
> c      APHI1=((1.D0-XK2*0.5D0)*ELK-ELE)/XKRHO12
> cc
> cc  2nd loop:
> 
> c         P=(RRC2+RHOS)**2+ZS**2+DD2**2
> c         XK2=4.D0*RRC2*RHOS/P
> c         XK=SQRT(XK2)
> c         XKRHO12=XK*SQRT(RHOS)     !    NB#4, P.3
> cC
> c      XK2S=1.D0-XK2
> c      DL=DLOG(1.D0/XK2S)
> c      ELK=1.38629436112d0+XK2S*(0.09666344259D0+XK2S*(0.03590092383+
> c     *     XK2S*(0.03742563713+XK2S*0.01451196212))) +DL*
> c     *     (0.5D0+XK2S*(0.12498593597D0+XK2S*(0.06880248576D0+
> c     *     XK2S*(0.03328355346D0+XK2S*0.00441787012D0))))
> c      ELE=1.D0+XK2S*(0.44325141463D0+XK2S*(0.0626060122D0+XK2S*
> c     *      (0.04757383546D0+XK2S*0.01736506451D0))) +DL*
> c     *     XK2S*(0.2499836831D0+XK2S*(0.09200180037D0+XK2S*
> c     *       (0.04069697526D0+XK2S*0.00526449639D0)))
> cC
> c      APHI2=((1.D0-XK2*0.5D0)*ELK-ELE)/XKRHO12
> 
> c      APPRC=A1*APHI1+A2*APHI2
> c      IF (PROX) APPRC=APPRC*SINT/SINT1   !   LINEAR INTERPOLATION, IF TOO CLOSE TO THE Z-AXIS
> cC
> c      RETURN
> c      END
> cC
2192c2195
<          SUBROUTINE PRC_QUAD (X,Y,Z,BX,BY,BZ)
---
>          SUBROUTINE T01PRC_QUAD (X,Y,Z,BX,BY,BZ)
2248,2378c2251,2381
< c
< c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
< C
<       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/
< 
<         SINT2=SINT**2
<         COST2=COST**2
<         SC=SINT*COST
<         ALPHA=SINT2/R
<         GAMMA=COST/R**2
< 
<         CALL FFS(ALPHA,AL1,DAL1,F,FA,FS)
<         D1=SC*F**XK1/((R/B1)**BE1+1.D0)
<         D2=D1*COST2
< 
<         CALL FFS(ALPHA,AL2,DAL2,F,FA,FS)
<         D3=SC*FS**XK2/((R/B2)**BE2+1.D0)
<         D4=D3*COST2
< 
<         CALL FFS(ALPHA,AL3,DAL3,F,FA,FS)
<         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
< 
<         CALL FFS(ALPHA,AL6,DAL6,F,FA,FS)
<         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
<         RETURN
<         END
< c
< C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
< C
<         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/
< 
<         SINT2=SINT**2
<         COST2=COST**2
<         SC=SINT*COST
<         ALPHA=SINT2/R
<         GAMMA=COST/R**2
< 
<         CALL FFS(ALPHA,AL1,DAL1,F,FA,FS)
<         D1=F**XK1/((R/B1)**BE1+1.D0)
<         D2=D1*COST2
< 
<         CALL FFS(ALPHA,AL2,DAL2,F,FA,FS)
<         D3=FA**XK2/R**BE2
<         D4=D3*COST2
< 
<         CALL FFS(ALPHA,AL3,DAL3,F,FA,FS)
<         D5=FS**XK3*ALPHA**XK4/((R/B3)**BE3+1.D0)
<         D6=D5*COST2
< 
<         CALL FFS(GAMMA,0.D0,DG1,F,FA,FS)
<         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
<        RETURN
<        END
---
> cc
> cc&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
> cC
> c      DOUBLE PRECISION FUNCTION BR_PRC_Q (R,SINT,COST)
> cC
> cCalculates the radial component of the "quadrupole" part of the model partial ring current.
> cC
> c      IMPLICIT  REAL * 8  (A - H, O - Z)
> 
> c      DATA A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11,A12,A13,A14,A15,A16,A17,   ! ALL LINEAR PARAMETERS HERE
> c     * A18,XK1,AL1,DAL1,B1,BE1,XK2,AL2,DAL2,B2,BE2,XK3,XK4,AL3,DAL3,B3,  ! WERE MULTIPLIED BY 0.1,
> c     * BE3,AL4,DAL4,DG1,AL5,DAL5,DG2,C1,C2,C3,AL6,DAL6,DRM/-21.2666329,  ! SO THAT THEY CORRESPOND TO P_0=1 nPa,
> c     *32.24527521,-6.062894078,7.515660734,233.7341288,-227.1195714,     ! RATHER THAN THE ORIGINAL VALUE OF 10 nPa
> c     *8.483233889,16.80642754,-24.63534184,9.067120578,-1.052686913,     ! ASSUMED IN THE BIOT-SAVART INTEGRAL
> c     *-12.08384538,18.61969572,-12.71686069,47017.35679,-50646.71204,
> c     *7746.058231,1.531069371,2.318824273,.1417519429,.6388013110E-02,
> c     *5.303934488,4.213397467,.7955534018,.1401142771,.2306094179E-01,
> c     *3.462235072,2.568743010,3.477425908,1.922155110,.1485233485,
> c     *.2319676273E-01,7.830223587,8.492933868,.1295221828,.01753008801,
> c     *.01125504083,.1811846095,.04841237481,.01981805097,6.557801891,
> c     *6.348576071,5.744436687,.2265212965,.1301957209,.5654023158/
> 
> c        SINT2=SINT**2
> c        COST2=COST**2
> c        SC=SINT*COST
> c        ALPHA=SINT2/R
> c        GAMMA=COST/R**2
> 
> c        CALL T01FFS(ALPHA,AL1,DAL1,F,FA,FS)
> c        D1=SC*F**XK1/((R/B1)**BE1+1.D0)
> c        D2=D1*COST2
> 
> c        CALL T01FFS(ALPHA,AL2,DAL2,F,FA,FS)
> c        D3=SC*FS**XK2/((R/B2)**BE2+1.D0)
> c        D4=D3*COST2
> 
> c        CALL T01FFS(ALPHA,AL3,DAL3,F,FA,FS)
> c        D5=SC*(ALPHA**XK3)*(FS**XK4)/((R/B3)**BE3+1.D0)
> c        D6=D5*COST2
> 
> c        ARGA=((ALPHA-AL4)/DAL4)**2+1.D0
> c        ARGG=1.D0+(GAMMA/DG1)**2
> 
> c        D7=SC/ARGA/ARGG
> c        D8=D7/ARGA
> c        D9=D8/ARGA
> c        D10=D9/ARGA
> 
> c        ARGA=((ALPHA-AL5)/DAL5)**2+1.D0
> c        ARGG=1.D0+(GAMMA/DG2)**2
> 
> c        D11=SC/ARGA/ARGG
> c        D12=D11/ARGA
> c        D13=D12/ARGA
> c        D14=D13/ARGA
> 
> 
> c        D15=SC/(R**4+C1**4)
> c        D16=SC/(R**4+C2**4)*COST2
> c        D17=SC/(R**4+C3**4)*COST2**2
> 
> c        CALL T01FFS(ALPHA,AL6,DAL6,F,FA,FS)
> c        D18=SC*FS/(1.D0+((R-1.2D0)/DRM)**2)
> 
> c        BR_PRC_Q=A1*D1+A2*D2+A3*D3+A4*D4+A5*D5+A6*D6+A7*D7+A8*D8+A9*D9+
> c     *  A10*D10+A11*D11+A12*D12+A13*D13+A14*D14+A15*D15+A16*D16+A17*D17+
> c     *   A18*D18
> cC
> c        RETURN
> c        END
> cc
> cC%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> cC
> c        DOUBLE PRECISION FUNCTION BT_PRC_Q (R,SINT,COST)
> cC
> cCalculates the Theta component of the "quadrupole" part of the model partial ring current.
> cC
> c        IMPLICIT  REAL * 8  (A - H, O - Z)
> 
> c      DATA A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11,A12,A13,A14,A15,A16,A17,  ! ALL LINEAR PARAMETERS HERE
> c     *XK1,AL1,DAL1,B1,BE1,XK2,AL2,DAL2,BE2,XK3,XK4,AL3,DAL3,B3,BE3,AL4, ! WERE MULTIPLIED BY 0.1,
> c     *DAL4,DG1,AL5,DAL5,DG2,C1,C2,C3/12.74640393,-7.516393516,          ! SO THAT THEY CORRESPOND TO P_0=1 nPa,
> c     *-5.476233865,3.212704645,-59.10926169,46.62198189,-.01644280062,  ! RATHER THAN THE ORIGINAL VALUE OF 10 nPa
> c     *.1234229112,-.08579198697,.01321366966,.8970494003,9.136186247,   ! ASSUMED IN THE BIOT-SAVART INTEGRAL
> c     *-38.19301215,21.73775846,-410.0783424,-69.90832690,-848.8543440,
> c     *1.243288286,.2071721360,.05030555417,7.471332374,3.180533613,
> c     *1.376743507,.1568504222,.02092910682,1.985148197,.3157139940,
> c     *1.056309517,.1701395257,.1019870070,6.293740981,5.671824276,
> c     *.1280772299,.02189060799,.01040696080,.1648265607,.04701592613,
> c     *.01526400086,12.88384229,3.361775101,23.44173897/
> 
> c        SINT2=SINT**2
> c        COST2=COST**2
> c        SC=SINT*COST
> c        ALPHA=SINT2/R
> c        GAMMA=COST/R**2
> 
> c        CALL T01FFS(ALPHA,AL1,DAL1,F,FA,FS)
> c        D1=F**XK1/((R/B1)**BE1+1.D0)
> c        D2=D1*COST2
> 
> c        CALL T01FFS(ALPHA,AL2,DAL2,F,FA,FS)
> c        D3=FA**XK2/R**BE2
> c        D4=D3*COST2
> 
> c        CALL T01FFS(ALPHA,AL3,DAL3,F,FA,FS)
> c        D5=FS**XK3*ALPHA**XK4/((R/B3)**BE3+1.D0)
> c        D6=D5*COST2
> 
> c        CALL T01FFS(GAMMA,0.D0,DG1,F,FA,FS)
> c        FCC=(1.D0+((ALPHA-AL4)/DAL4)**2)
> c        D7 =1.D0/FCC*FS
> c        D8 =D7/FCC
> c        D9 =D8/FCC
> c        D10=D9/FCC
> 
> c        ARG=1.D0+((ALPHA-AL5)/DAL5)**2
> c        D11=1.D0/ARG/(1.D0+(GAMMA/DG2)**2)
> c        D12=D11/ARG
> c        D13=D12/ARG
> c        D14=D13/ARG
> 
> c        D15=1.D0/(R**4+C1**2)
> c        D16=COST2/(R**4+C2**2)
> c        D17=COST2**2/(R**4+C3**2)
> cC
> c        BT_PRC_Q=A1*D1+A2*D2+A3*D3+A4*D4+A5*D5+A6*D6+A7*D7+A8*D8+A9*D9+
> c     *   A10*D10+A11*D11+A12*D12+A13*D13+A14*D14+A15*D15+A16*D16+A17*D17
> cC
> c       RETURN
> c       END
2382c2385
<        SUBROUTINE FFS(A,A0,DA,F,FA,FS)
---
>        SUBROUTINE T01FFS(A,A0,DA,F,FA,FS)
2395c2398
<          SUBROUTINE RC_SHIELD (A,PS,X_SC,X,Y,Z,BX,BY,BZ)
---
>          SUBROUTINE T01RC_SHIELD (A,PS,X_SC,X,Y,Z,BX,BY,BZ)
2542c2545
<       SUBROUTINE DIPOLE01 (PS,X,Y,Z,BX,BY,BZ)
---
>       SUBROUTINE T01DIPOLE01 (PS,X,Y,Z,BX,BY,BZ)
