c                                 PERMAG
c
c                          Version of Apr 16, 1987
c                          Written by Z.Q. Xie and T.A. Antaya
c
c                          NSCL ,  MSU
c                          E. Lansing, MI  48824-1321
c
c
c-------------------------------------------------------------------------------
c-----------
c
c      PERMAG will calculate the B field components of dipole, quadrupole,
c  hexapole, octupole and higher order multipoles that are made of geometrically
c  symmetric hexagonal or trapezoidal or rectangular permanent magnet blocks.
c  The main assumption made is that there is a high anisotropy in the
c  magnetization direction. PERMAG does also the case of the field produced by
c  a single permanent magnet. PERMAG calculates the multipole field about a
c  round aperture, but with some straightforward modifications, it could also
c  calculate rectangular or irregular multipole aperture fields as well as
c  calculations of multipole field which the multipole is made of combination
c  of rectangular and trapezoidal magnet blocks.
c
c-------------------------------------------------------------------------------
c
c                            DATA INPUT
c*******************************************************************************
c
c     PERMAG runs interactively. SUMMARY of all input prompts, as well as a
c  set of examples (test) data (all the geomtry dimensions have the same unit)
c  are given bellow:
c
c SUMMARY:
c
c  PROMPT 1  ' Type  N and M '
c
c           2N is the multipole order number and M is number of pieces. For
c           example: N=3, M=12 is a 12 piece hexapole. If M=1 , then a sigle
c           permanent magnet calculation is assumed  and 2N will be set to be 1.
c
c
c  PROMPT 2  ' Type R, A, B, D, E , F , C and AN1 '
c
c           R is the aperture radius, say , the distance from the aperture
c           center to a pole tip.  A, B, D , E , F and C are the physical
c           dimensions of the a single piece permanent magnet ( a multipole
c           consists of M pieces such magnets ). AN1 is the angle between the
c           magnetization of M=1 piece and the X-axis which bisects the M=1
c           piece.  This geometry is given in the report " THE PERMAG -- FOR
c           THE CALCULATION OF THE AIR CORE 3D MULTIPOLE FIELD PRODUCED BY
c           ORIENTED PERMANENT MAGNETS ".
c
c
c  PROMPT 3  ' Type  Brms (in kG) in the order of M = 1 ttroug M '
c
c           The M magnets could have M different residual inductions.
c           If M magnets have the same residual inductions, input M*Brm
c
c
c  PROMPT 4  ' Type  IC '
c
c           This selects the type of calculation to be performed. PERMAG has
c           5 different choices:
c
c           IC = 1 , Field calculations at constant Z. Specifically, the code
c                    calculates Br , Bphi , Bz , and the total strength Bt ,
c                    all vs r and phi. Z = 0 is assumed to be the multipole
c                    midplane.
c
c                  PROMPT 5  ' Type Z '
c
c                  PROMPT 6  ' Type  r_min , r_step , number of steps '
c                           r_min is the minimum r in which the calculation
c                           will start with in the r dimension,  r_step is the
c                           step size in the r dimension. Number of steps are
c                           the steps which will be stepped in the r dimension.
c
c                  PROMPT 7  ' Type Phi_min , Phi_step , number of steps '
c                           phi = 0 corresponds to the middle of M = 1 piece.
c                           Phi_min is the minimum phi in which the calculation
c                           will start with in the phi dimension. Phi_step is
c                           the step size in the phi dimension, and number of
c                           steps are the steps which will be stepped in the
c                           phi dimension.
c
c
c           IC = 2 , Field calculations at fixed azmuth phi = constant.
c                    Specifically, the code calculates Br , Bphi , Bz
c                    and the total strength Bt , all vs r and Z .
c
c                  PROMPT 5  ' Type phi (in degrees) '
c
c                  PROMPT 6  ' Type r_min , r_step , number of steps '
c
c                  PROMPT 7  ' Type Z_min , Z_step , number of steps '
c                           Z_min is the minimum Z in which the calculation
c                           will start with in the Z dimension. Z_step is the
c                           step size in the Z dimension and number of steps
c                           are the steps which will be stepped in the Z
c                           dimension.
c
c
c           IC = 3 , Field calculations at fixed r. Specifically, the code
c                    calculates Br , Bphi , Bz and the total strength Bt ,
c                    all vs phi and Z.
c
c                  PROMPT 5  ' Type r '
c
c                  PROMPT 6  ' Type Z_min , Z_step , number of steps '
c                           See definitions above.
c
c                  PROMPT 7  ' Type Phi_min , Phi_step , number of steps '
c                           See definitions above.
c
c
c           IC = 4 , The code calculates the total field strength vs r at a
c                    constant phi and constant Z.
c
c                  PROMPT 5  ' Type  Z ,r_min , r_step ,number of steps in r '
c
c                  PROMPT 6  ' Type  phi ( in degrees ) '
c
c
c
c           IC = 5 , Field calculation at a single point (r,phi,Z), and
c                    Br , Bphi , Bz and Bt are returned to the screen.
c
c                  PROMPT 5  ' Type r , phi , Z  '
c
c
c
c         ---------------------------------------------------------------------
c
c EXAMPLE:
c
c       A hexapole ( N = 3 ) made of 12 ( M = 12 ) trapezoidal pieces, all
c       having the same residual induction Brm = 8.7 kG, the report on this
c       code includes the output files for this example calculation in
c       Appendix C.
c
c1     3,12                                   !  N , M
c2     5,4,2,3.5,0,0,50,0                     !  R , A , B , D , E , F , C , AN1
c3     12*8.7   (or 8.7,8.7,8.7,8.7,........) !  Brm(1), Brm(2)........Br(12)
c4     4                                      !  IC
c5     0,0,.25,20                             !  Z,r_min , r_step, NS
c6     30. ( or any phi )                     !  Phi ( in degrees )
c
c*******************************************************************************
c
c       SUMMARY OF DATA OUTPUT
c*******************************************************************************
c
c      The output file named " PERMAG.OUT " , in which all the input dimensions
c      informations are recorded, such as N , M , R , A , B , D , E , F ,C ,
c      AN1 ,Brms and the fixed variables. The strength B is in kG unit and all
c      the physical dimensions have the same units.
c
c*******************************************************************************
c*******************************************************************************
*******
c
c   Special Note for M = 1 :
c
c       When M =1 , a single permanent magnet only , one could, if one wants to,
c     translate the origin of the current coordinate to the center of the
c     single piece by simply inputting R = -A/2.
c
c*******************************************************************************
c******
c
c
      character*1  REPLY,CHANGE
      character*7  residual
      real*8 A,B,D,C,E,F,R,AN,r_min,r_step,phi_min,SI34,CO34,SI56,CO56,Z
      real*8 phi_step,Z_min,Z_step,X(500),Y(500)
      real*8 PHID,PHIM,PHIM1,PHI,theta34,theta56,easy_axis,P
      real*8 HX(100),HY(100),HZ(100),BR(50)
      real*8 OX1,OX2,OX3,OX4,OX5,OX6,OX7
      real*8 OX8,OX9,OX10,OX11,OX12,OX13,OX14,OX15,OX16
      real*8 OX17,OX18,OX19,OX20,OX21,OX22,OX23,OX24
      real*8 OY1,OY2,OY3,OY4,OY5,OY6,OY7,OY8
      real*8 OY9,OY10,OY11,OY12,OY13,OY14,OY15,OY16
      real*8 OP1,OP2,OP3,OP4,OP5,OP6,OP7,OP8,OP9,OP10,OP11,OP12,OP13
      real*8 OP14,OP15,OP16,OP17,OP18,OP19,OP20,OP21,OP22,OP23,OP24
      real*8 TP1,TP2,TP3,TP4,TP5,TP6,TP7,TP8,TP9,TP10,TP11,TP12,TP13
      real*8 TP14,TP15,TP16
      real*8 BTY,BTZ,BTX3,BTX4,BTX5,BTX6,BTZ3,BTZ4,BTZ5,BTZ6
      real*8 BY1,BY2,BY3,BY4,BY5,BY6,BY7,BY8
      real*8 BZ1,BZ2,BZ3,BZ4,BZ5,BZ6,BZ7,BZ8
      real*8 BYY1,BYY2,BYY3,BYY4,BYY5,BYY6,BYY7,BYY8
      real*8 BZZ1,BZZ2,BZZ3,BZZ4,BZZ5,BZZ6,BZZ7,BZZ8
      real*8 BX1,BX2,BX3,BX4,BX5,BX6,BX7,BX8,BX9,BX10,BX11,BX12
      real*8 BXX1,BXX2,BXX3,BXX4,BXX5,BXX6,BXX7,BXX8,BXX9,BXX10,BXX11
      real*8 BXX12,BLZ1,BLZ2,BLZ3,BLZ4,BLZ5,BLZ6,BLZ7,BLZ8,BLZ9,BLZ10
      real*8 BLZ11,BLZ12,BLZ13,BLZ14,BLZ15,BLZ16,BL14,BL15,BL16
      real*8 BL1,BL2,BL3,BL4,BL5,BL6,BL7,BL8,BL9,BL10,BL11,BL12,BL13
      real*8 I1,I2,I3,Z0,BZL1(100),BZL2(100),BZL3(100)
      real*8 BRAL,HHZL,BRAT,HHZT,BTLH,BPERP
cc
c  X,Y,Z contributions from surface 1
cc
      OX1(X,Y,Z)=atan((Y-B/2.)*(Z-C/2.)/((X-R)*((X-R)**2.+
     1(Y-B/2.)**2+(Z-C/2.)**2.)**.5))
      OX2(X,Y,Z)=atan(-(Y-B/2.)*(Z+C/2.)/((X-R)*((X-R)**2.+
     1(Y-B/2.)**2+(Z+C/2.)**2.)**.5))
      OX3(X,Y,Z)=atan((Y+B/2.)*(Z+C/2.)/((X-R)*((X-R)**2.+
     1(Y+B/2.)**2.+(Z+C/2.)**2.)**.5))
      OX4(X,Y,Z)=atan(-(Y+B/2.)*(Z-C/2.)/((X-R)*((X-R)**2.+
     1(Y+B/2.)**2.+(Z-C/2.)**2)**.5))
 
      BY1(X,Y,Z)=(Z+C/2.)+((Z+C/2.)**2.+(X-R)**2.+(Y-B/2.)**2.)**.5
      BY2(X,Y,Z)=(Z-C/2.)+((Z-C/2.)**2.+(X-R)**2.+(Y+B/2.)**2.)**.5
      BY3(X,Y,Z)=(Z-C/2.)+((Z-C/2.)**2.+(X-R)**2.+(Y-B/2.)**2.)**.5
      BY4(X,Y,Z)=(Z+C/2.)+((Z+C/2.)**2.+(X-R)**2.+(Y+B/2.)**2.)**.5
      BZ1(X,Y,Z)=(Y+B/2.)+((Z-C/2.)**2.+(X-R)**2.+(Y+B/2.)**2.)**.5
      BZ2(X,Y,Z)=(Y-B/2.)+((Z+C/2.)**2.+(X-R)**2.+(Y-B/2.)**2.)**.5
      BZ3(X,Y,Z)=(Y-B/2.)+((Z-C/2.)**2.+(X-R)**2.+(Y-B/2.)**2.)**.5
      BZ4(X,Y,Z)=(Y+B/2.)+((Z+C/2.)**2.+(X-R)**2.+(Y+B/2.)**2.)**.5
cc
c
c  X,Y,Z contributions from surface 2
c
      OX5(X,Y,Z)=atan((Y-F/2.)*(Z-C/2.)/((X-R-A-E)*((X-R-A-E)*
     1*2.+(Y-F/2.)**2.+(Z-C/2.)**2)**.5))
      OX6(X,Y,Z)=atan(-(Y-F/2.)*(Z+C/2.)/((X-R-A-E)*((X-R-A-E)*
     1*2.+(Y-F/2.)**2.+(Z+C/2.)**2.)**.5))
      OX7(X,Y,Z)=atan((Y+F/2.)*(Z+C/2.)/((X-R-A-E)*((X-R-A-E)*
     1*2.+(Y+F/2.)**2.+(Z+C/2.)**2.)**.5))
      OX8(X,Y,Z)=atan(-(Y+F/2.)*(Z-C/2.)/((X-R-A-E)*((X-R-A-E)
     1**2.+(Y+F/2.)**2.+(Z-C/2.)**2.)**.5))
      BY5(X,Y,Z)=(Z+C/2.)+((Z+C/2.)**2.+(X-R-A-E)**2.+(Y+F/2.)**2.)**.5
      BY6(X,Y,Z)=(Z-C/2.)+((Z-C/2.)**2.+(X-R-A-E)**2.+(Y-F/2.)**2.)**.5
      BY7(X,Y,Z)=(Z+C/2.)+((Z+C/2.)**2.+(Y-F/2.)**2.+(X-R-A-E)**2.)**.5
      BY8(X,Y,Z)=(Z-C/2.)+((Z-C/2.)**2.+(X-R-A-E)**2.+(Y+F/2.)**2.)**.5
      BZ5(X,Y,Z)=(Y-F/2.)+((Z-C/2.)**2.+(X-A-R-E)**2.+(Y-F/2.)**2.)**.5
      BZ6(X,Y,Z)=(Y+F/2.)+((Z+C/2.)**2.+(X-R-A-E)**2.+(Y+F/2.)**2.)**.5
      BZ7(X,Y,Z)=(Y+F/2.)+((Z-C/2.)**2.+(Y+F/2.)**2.+(X-R-A-E)**2.)**.5
      BZ8(X,Y,Z)=(Y-F/2.)+((Z+C/2.)**2.+(X-R-A-E)**2.+(Y-F/2.)**2.)**.5
cc
c
c  X,Y,Z contributions from surface 3
cc
      OX9(X,Y,Z)=atan((Z-C/2.)*((X-R)*CO34-(Y+B/2.)*SI34)/(((X-R)*
     1SI34+(Y+B/2.)*CO34)*((X-R)**2.+(Y+B/2.)**2.+(Z-C/2.)**2.)**.5))
      OX10(X,Y,Z)=atan(-(Z-C/2.)*((X-R-A)*CO34-(Y+D/2.)*SI34)/(((X-R
     1-A)*SI34+(Y+D/2.)*CO34)*((X-R-A)**2.+(Y+D/2.)**2.+(Z-C/2.)**2.
     1)**.5))
      OX11(X,Y,Z)=atan((Z+C/2.)*((X-R-A)*CO34-(Y+D/2.)*SI34)/(((X-R-A)*
     1SI34+(Y+D/2.)*CO34)*((X-R-A)**2.+(Y+D/2.)**2.+(Z+C/2.)**2.)**.5))
      OX12(X,Y,Z)=atan(-(Z+C/2.)*((X-R)*CO34-(Y+B/2.)*SI34)/(((X-R)*
     1SI34+(Y+B/2.)*CO34)*((X-R)**2.+(Y+B/2.)**2.+(Z+C/2.)**2.)**.5))
      BX1(X,Y,Z)=(Z+C/2.)+((Z+C/2.)**2.+(Y+D/2.)**2.+(X-R-A)**2.)**.5
      BX2(X,Y,Z)=(Z-C/2.)+((Z-C/2.)**2.+(X-R)**2.+(Y+B/2.)**2.)**.5
      BX3(X,Y,Z)=(Z+C/2.)+((Z+C/2.)**2.+(X-R)**2.+(Y+B/2.)**2.)**.5
      BX4(X,Y,Z)=(Z-C/2.)+((Z-C/2.)**2.+(X-R-A)**2.+(Y+D/2.)**2.)**.5
 
      OY1(X,Y,Z)=atan(-(Z-C/2.)*((X-R)*CO34-(Y+B/2.)*SI34)/(((X-R)*
     1SI34+(Y+B/2.)*CO34)*((X-R)**2.+(Y+B/2.)**2.+(Z-C/2.)**2.)**.5))
      OY2(X,Y,Z)=atan((Z-C/2.)*((X-R-A)*CO34-(Y+D/2.)*SI34)/(((X-R-A)*
     1SI34+(Y+D/2.)*CO34)*((X-R-A)**2.+(Y+D/2.)**2.+(Z-C/2.)**2.)**.5))
      OY3(X,Y,Z)=atan(-(Z+C/2.)*((X-R-A)*CO34-(Y+D/2.)*SI34)/(((X-R-A)*
     1SI34+(Y+D/2.)*CO34)*((X-R-A)**2.+(Y+D/2.)**2.+(Z+C/2.)**2.)**.5))
      OY4(X,Y,Z)=atan((Z+C/2.)*((X-R)*CO34-(Y+B/2.)*SI34)/(((X-R)*
     1SI34+(Y+B/2.)*CO34)*((X-R)**2.+(Y+B/2.)**2.+(Z+C/2.)**2.)**.5))
      BLZ1(X,Y,Z)=((Y+B/2.)*SI34-(X-R)*CO34)+((Z+C/2.)**2.+(Y+B/2.)*
     1*2.+(X-R)**2.)**.5
      BLZ2(X,Y,Z)=((Y+D/2.)*SI34-(X-R-A)*CO34)+((Z-C/2.)**2.+(X-R
     1-A)**2.+(Y+D/2.)**2.)**.5
      BLZ3(X,Y,Z)=((Y+B/2.)*SI34-(X-R)*CO34)+((Z-C/2.)
     1**2.+(X-R)**2.+(Y+B/2.)**2.)**.5
      BLZ4(X,Y,Z)=((Y+D/2.)*SI34-(X-R-A)*CO34)+
     1((Z+C/2.)**2.+(X-A-R)**2.+(Y+D/2.)**2.)**.5
cc
c
c  X,Y,Z contributions from surface 4
cc
      OX13(X,Y,Z)=atan((Z-C/2.)*((X-R)*CO34+(Y-B/2.)*SI34)/(((X-R)*SI34
     1-(Y-B/2.)*CO34)*((X-R)**2.+(Y-B/2.)**2.+(Z-C/2.)**2.)**.5))
      OX14(X,Y,Z)=atan(-(Z-C/2.)*((X-R-A)*CO34+(Y-D/2.)*SI34)/(((X-R-A)
     1*SI34-(Y-D/2.)*CO34)*((X-R-A)**2.+(Y-D/2.)**2.+(Z-C/2.)**2.)**.5))
      OX15(X,Y,Z)=atan((Z+C/2.)*((X-R-A)*CO34+(Y-D/2.)*SI34)/(((X-R-A)*
     1SI34-(Y-D/2.)*CO34)*((X-R-A)**2.+(Y-D/2.)**2.+(Z+C/2.)**2.)**.5))
      OX16(X,Y,Z)=atan(-(Z+C/2.)*((X-R)*CO34+(Y-B/2.)*SI34)/(((X-R)*
     1SI34-(Y-B/2.)*CO34)*((X-R)**2.+(Y-B/2.)**2.+(Z+C/2.)**2.)**.5))
 
      OY5(X,Y,Z)=atan((Z-C/2.)*((X-R)*CO34+(Y-B/2.)*SI34)/(((X-R)*SI34-
     1(Y-B/2.)*CO34)*((X-R)**2.+(Y-B/2.)**2.+(Z-C/2.)**2.)**.5))
      OY6(X,Y,Z)=atan(-(Z-C/2.)*((X-R-A)*CO34+(Y-D/2.)*SI34)/(((X-R-A)*
     1SI34-(Y-D/2.)*CO34)*((X-R-A)**2.+(Y-D/2.)**2.+(Z-C/2.)**2.)**.5))
      OY7(X,Y,Z)=atan((Z+C/2.)*((X-R-A)*CO34+(Y-D/2.)*SI34)/(((X-R-A)*
     1SI34-(Y-D/2.)*CO34)*((X-R-A)**2.+(Y-D/2.)**2.+(Z+C/2.)**2.)**.5))
      OY8(X,Y,Z)=atan(-(Z+C/2.)*((X-R)*CO34+(Y-B/2.)*SI34)/(((X-R)*
     1SI34-(Y-B/2.)*CO34)*((X-R)**2.+(Y-B/2.)**2.+(Z+C/2.)**2.)**.5))
 
      BX5(X,Y,Z)=(Z+C/2.)+((Z+C/2.)**2.+(Y-D/2.)**2.+(X-R-A)**2.)**.5
      BX6(X,Y,Z)=(Z-C/2.)+((Z-C/2.)**2.+(X-R)**2.+(Y-B/2.)**2.)**.5
      BX7(X,Y,Z)=(Z+C/2.)+((Z+C/2.)**2.+(X-R)**2.+(Y-B/2.)**2.)**.5
      BX8(X,Y,Z)=(Z-C/2.)+((Z-C/2.)**2.+(X-R-A)**2.+(Y-D/2.)**2.)**.5
 
      BLZ5(X,Y,Z)=-(Y-D/2.)*SI34-(X-R-A)*CO34+((Z+C/2.)**2.+(X-A
     1-R)**2.+(Y-D/2.)**2.)**.5
      BLZ6(X,Y,Z)=-(Y-B/2.)*SI34-(X-R)*CO34+((Z-C/2.)**2.+(X-R)*
     1*2.+(Y-B/2.)**2.)**.5
      BLZ7(X,Y,Z)=-(Y-B/2.)*SI34-(X-R)*CO34+((Z+C/2.)**2.
     1+(X-R)**2.+(Y-B/2.)**2.)**.5
      BLZ8(X,Y,Z)=-(Y-D/2.)*SI34-(X-R-A)*CO34+((Z-
     1C/2.)**2.+(X-A-R)**2.+(Y-D/2.)**2.)**.5
CC
c  X,Y,Z contributions from surface 5
c  There are some terms are the same as those in the surface 3, so they are
c   omitted here.
cc
      OX17(X,Y,Z)=atan((Z-C/2.)*((X-R-A)*CO56-(Y+D/2.)*SI56)/(
     1((X-R-A)*SI56+(Y+D/2.)*CO56)*((X-R-A)**2.+(Y+D/2.)**2.+
     1(Z-C/2.)**2.)**.5))
      OX18(X,Y,Z)=atan(-(Z-C/2.)*((X-R-A-E)*CO56-(Y+F/2.)*SI56)
     1/(((X-R-A-E)*SI56+(Y+F/2.)*CO56)*((X-R-A-E)**2.+(Y+F/2.)*
     1*2.+(Z-C/2.)**2.)**.5))
      OX19(X,Y,Z)=atan((Z+C/2.)*((X-R-A-E)*CO56-(Y+F/2.)*SI56)/
     1(((X-R-A-E)*SI56+(Y+F/2.)*CO56)*((X-R-A-E)**2.+(Y+F/2.)*
     1*2.+(Z+C/2.)**2.)**.5))
      OX20(X,Y,Z)=atan(-(Z+C/2.)*((X-R-A)*CO56-(Y+D/2.)*SI56)/
     1(((X-R-A)*SI56+(Y+D/2.)*CO56)*((X-R-A)**2.+(Y+D/2.)**2.+
     1(Z+C/2.)**2.)**.5))
      BX9(X,Y,Z)=(Z-C/2.)+((Z-C/2.)**2.+(X-R-A-E)**2.+(Y+F/2.
     1)**2.)**.5
      BX10(X,Y,Z)=(Z+C/2.)+((Z+C/2.)**2.+(X-R-A-E)**2.+(Y+F/2.
     1)**2.)**.5
 
      OY9(X,Y,Z)=atan(-(Z-C/2.)*((X-R-A)*CO56-(Y+D/2.)*SI56)/(
     1((X-R-A)*SI56+(Y+D/2.)*CO56)*((X-R-A)**2.+(Y+D/2.)**2.+
     1(Z-C/2.)**2.)**.5))
      OY10(X,Y,Z)=atan((Z-C/2.)*((X-R-A-E)*CO56-(Y+F/2.)*SI56)
     1/(((X-R-A-E)*SI56+(Y+F/2.)*CO56)*((X-R-A-E)**2.+(Y+F/2.)*
     1*2.+(Z-C/2.)**2.)**.5))
      OY11(X,Y,Z)=atan(-(Z+C/2.)*((X-R-A-E)*CO56-(Y+F/2.)*SI56
     1)/(((X-R-A-E)*SI56+(Y+F/2.)*CO56)*((X-R-A-E)**2.+(Y+F/2.
     1)**2.+(Z+C/2.)**2.)**.5))
      OY12(X,Y,Z)=atan((Z+C/2.)*((X-R-A)*CO56-(Y+D/2.)*SI56)/(
     1((X-R-A)*SI56+(Y+D/2.)*CO56)*((X-R-A)**2.+(Y+D/2.)**2.+
     1(Z+C/2.)**2.)**.5))
      BLZ9(X,Y,Z)=((Y+D/2.)*SI56-(X-R-A)*CO56)+((Z+C/2.)**2.+
     1(Y+D/2.)**2.+(X-R-A)**2.)**.5
      BLZ10(X,Y,Z)=((Y+F/2.)*SI56-(X-R-A-E)*CO56)+((Z-C/2.)*
     1*2.+(X-R-A-E)**2.+(Y+F/2.)**2.)**.5
      BLZ11(X,Y,Z)=((Y+D/2.)*SI56-(X-R-A)*CO56)+((Z-C/2.)
     1**2.+(X-R-A)**2.+(Y+D/2.)**2.)**.5
      BLZ12(X,Y,Z)=((Y+F/2.)*SI56-(X-R-A-E)*CO56)+
     1((Z+C/2.)**2.+(X-A-R-E)**2.+(Y+F/2.)**2.)**.5
cc
c
c  X,Y,Z contributions from surface 6
c  There are some terms are the same as those in the surface 4, so they are
c   omitted here.
cc
      OX21(X,Y,Z)=atan((Z-C/2.)*((X-R-A)*CO56+(Y-D/2.)*SI56)
     1/(((X-R-A)*SI56-(Y-D/2.)*CO56)*((X-R-A)**2.+(Y-D/2.)*
     1*2.+(Z-C/2.)**2.)**.5))
      OX22(X,Y,Z)=atan(-(Z-C/2.)*((X-R-A-E)*CO56+(Y-F/2.)*
     1SI56)/(((X-R-A-E)*SI56-(Y-F/2.)*CO56)*((X-R-A-E)**2.+
     1(Y-F/2.)**2.+(Z-C/2.)**2.)**.5))
      OX23(X,Y,Z)=atan((Z+C/2.)*((X-R-A-E)*CO56+(Y-F/2.)*
     1SI56)/(((X-R-A-E)*SI56-(Y-F/2.)*CO56)*((X-R-A-E)**2.
     1+(Y-F/2.)**2.+(Z+C/2.)**2.)**.5))
      OX24(X,Y,Z)=atan(-(Z+C/2.)*((X-R-A)*CO56+(Y-D/2.)*SI56
     1)/(((X-R-A)*SI56-(Y-D/2.)*CO56)*((X-R-A)**2.+(Y-D/2.)
     1**2.+(Z+C/2.)**2.)**.5))
 
      OY13(X,Y,Z)=atan((Z-C/2.)*((X-R-A)*CO56+(Y-D/2.)*SI56
     1)/(((X-R-A)*SI56-(Y-D/2.)*CO56)*((X-R-A)**2.+(Y-D/2.)
     1**2.+(Z-C/2.)**2.)**.5))
      OY14(X,Y,Z)=atan(-(Z-C/2.)*((X-R-A-E)*CO56+(Y-F/2.)*
     1SI56)/(((X-R-A-E)*SI56-(Y-F/2.)*CO56)*((X-R-A-E)**2.
     1+(Y-F/2.)**2.+(Z-C/2.)**2.)**.5))
      OY15(X,Y,Z)=atan((Z+C/2.)*((X-R-A-E)*CO56+(Y-F/2.)*
     1SI56)/(((X-R-A-E)*SI56-(Y-F/2.)*CO56)*((X-R-A-E)**2.
     1+(Y-F/2.)**2.+(Z+C/2.)**2.)**.5))
      OY16(X,Y,Z)=atan(-(Z+C/2.)*((X-R-A)*CO56+(Y-D/2.)*SI56
     1)/(((X-R-A)*SI56-(Y-D/2.)*CO56)*((X-R-A)**2.+(Y-D/2.)
     1**2.+(Z+C/2.)**2.)**.5))
 
      BX11(X,Y,Z)=(Z-C/2.)+((Z-C/2.)**2.+(X-R-A-E)**2.+
     1(Y-F/2.)**2.)**.5
      BX12(X,Y,Z)=(Z+C/2.)+((Z+C/2.)**2.+(Y-F/2.)**2.+
     1(X-R-A-E)**2.)**.5
 
      BLZ13(X,Y,Z)=-(Y-F/2.)*SI56-(X-R-A-E)*CO56+((Z+C/2.)*
     1*2.+(X-A-R-E)**2.+(Y-F/2.)**2.)**.5
      BLZ14(X,Y,Z)=-(Y-D/2.)*SI56-(X-R-A)*CO56+((Z-C/2.)**2.
     1+(X-R-A)**2.+(Y-D/2.)**2.)**.5
      BLZ15(X,Y,Z)=-(Y-D/2.)*SI56-(X-R-A)*CO56+((Z+C/2.)**2.
     1+(X-R-A)**2.+(Y-D/2.)**2.)**.5
      BLZ16(X,Y,Z)=-(Y-F/2.)*SI56-(X-R-A-E)*CO56+((Z-
     1C/2.)**2.+(X-A-R-E)**2.+(Y-F/2.)**2.)**.5
 
      P=3.141592653589793
CC
C   INPUT.
CC
cc
c  Read in B-field values for L-2 coils
cc
      READ(21,261)I1,I2,I3,Z0
      I1=I1/150.0
      I2=I2/150.0
      I3=I3/150.0
  261 FORMAT(F7.2,3(3X,F7.2))
      DO K=1,100
        READ(21,262,end=271)BZL1(K),BZL2(K),BZL3(K)
        BZL1(K)=BZL1(K)/1000.0
        BZL2(K)=BZL2(K)/1000.0
        BZL3(K)=BZL3(K)/1000.0
      END DO
  262 FORMAT(7X,3(3X,F7.2))
  271 CONTINUE
      open(10,file='PERMAG.OUT',status='new')
      write(10,99)
99    format(//,5x,' Magnet Specifications: ',/)
1     type *, ' Type    N  ,  M  '
      accept *,IN,M
      if((M.lt.IN).and.(M.ne.1)) then
60     write(6,61)M
61     format(5x,' Your input: M =',I4,' can not construct such',
     1 ' multipole.  Reenter N  and  M ')
       accept *,IN,M
      if((M.lt.IN).and.(M.ne.1)) goto 60
       end if
        N=2*IN
       write(10,11)IN,M
11    format(//,5x,' N = ',I2,10x,' M = ',I2)
2     type *, ' Type    R , A , B , D , E , F ,C , AN1 '
      accept *,R,A,B,D,E,F,C,AN1
      AN=AN1*P/180.
       write(10,22)R,A,B,D,E,F,C,AN1
22     format(/,' R =',f8.3,2x,' A =',f8.3,2x,' B =',f8.3,2x,
     1 ' D =',f8.3,2x,//,' E =',f8.3,2x,' F =',f8.3,2x,' C =',
     1 f10.2,2x,' AN1 =',f8.2,/)
3     type *, ' Type    Brms ( in kGs ) in the order of 1 to M '
      read (5,*)(BR(I),I=1,M)
      do k=1,M
         residual(1:7)='Brm(  )'
         write(residual(5:6),32)k
         write(10,33)residual,Br(k)
32       format(I2)
33       format(6x,A7,' =',f6.3)
      enddo
      if(E.eq.0.) then
            F=D
      endif
4     type *, ' Type    IC '
      accept *,IC
      if(IC.ge.6) then
        type *,'  IC could not be larger than 5, your IC is not valid.'
        type *,'     Reinput IC '
        accept *, IC
       endif
      write (10,88)IC
88    format(///,11x,'IC =',I2)
      NS=0
      NSM=0
      goto (111,112,113,114,115), IC
111      type *,' Type   Z '
         accept *, Z1
       write(10,311)Z1
311    format(/,10x,' Z =',f6.1,//)
       write(10,312)
312    format(/,9x,'r',8x,'phi',13x,'Br',9x,'Bphi',9x,'Bz',10x,'Bt',/)
         type *,' Type   r_min  , r_step, number of steps '
         accept *, r_min,r_step,NS
         NN=NS+1
         type *,' Type   phi_min , phi_step , number of steps '
         accept *, phi_min,phi_step,NSM
         MS=NSM+1
       goto 199
112      type *,' Type   phi ( in degrees ) '
         accept *, PHI2
       write(10,322)PHI2
322    format(/,10x,' PHI =',f6.1,//)
       write(10,323)
323    format(/,9x,'r',9x,'Z',14x,'Br',9x,'Bphi',9x,'Bz',10x,'Bt',/)
         type *,' Type   r_min  , r_step, number of steps '
         accept *, r_min,r_step,NS
         NN=NS+1
         type *,' Type   Z_min , Z_step , number of steps '
         accept *, Z_min,Z_step,NSM
         MS=NSM+1
         goto 199
113      type *,' Type    r '
         Accept *, r3
       write(10,333)r3
333    format(/,10x,' r =',f6.1,//)
       write(10,334)
334    format(/,9x,'Z',8x,'phi',13x,'Br',9x,'Bphi',9x,'Bz',10x,'Bt',/)
         type *,' Type   Z_min , Z_step , number of steps '
         accept *, Z_min,Z_step,NS
         NN=NS+1
         type *,' Type   phi_min , phi_step , number of steps '
         accept *, phi_min,phi_step,NSM
         MS=NSM+1
           goto 199
114      type *,' Type  Z , r_min , r_step , number of steps in r '
         accept *, Z4,r_min,r_step,NS
         NN=NS+1
         MS=1
         type *,'  Type    phi ( in degrees )'
         accept *, PHI4
       write(10,343)Z4,PHI4
343    format(/,10x,' Z =',f6.1,10x,' PHI =',f6.1,//)
       write(10,342)
342    format(/,9x,'r',13x,'Br',9x,'Bphi',9x,'Bz',10x,'Bt',/)
         goto 199
115      type *,' Type    r , phi , Z '
         accept *,r5,PHI5,Z5
         write(10,355)
355    format(/,5x,'r',5x,'phi',5x,'Z',10x,'Br',10x,'Bphi',10x,
     1 'Bz',10x,'Bt',/)
         NN=1
         MS=1
cc
c   check the input to see if it is O.K. to construct the multipole.
cc
199   WRITE(10,201)NS,NSM
201   FORMAT(1X,2I3)
      OK=atan(B/2/R)-P/M
      KO=atan(D/2/(R+A))-P/M
      KOP=atan(F/2/(R+A+E))-P/M
      if((OK.gt.0.).or.(KO.gt.0.).or.(KOP.gt.0.)) then
      write (6,310)
310   format('  Dear Sir;'/
     1/'  Sorry, Your input dimensions are not right to construct'
     1//'        the multipole you want, try again please.' )
         goto 1
      endif
cc
c    " theta34 " is the side angle of the trapezoidal piece,
c    " theta56 " is the second side angle of the semmetrical polygon detail
c  refer to the report " AN ANALYTICAL CALCULATION OF MULTIPOLE FIELD PRODUCED
c  BY ORIENTED PERMANENT MAGNETS ".
c
200   theta34=atan((D-B)/2./A)
      if(E.ne.0.) then
        theta56=atan((F-D)/2./E)
        SI56=SIN(theta56)
        CO56=COS(theta56)
      endif
      SI34=SIN(theta34)
      CO34=COS(theta34)
      if (IC.eq.1) Z=Z1
      if (IC.eq.4) Z=Z4
      if(IC.eq.2) PHI=PHI2*P/180.
      if (IC.eq.3) Xr=r3
 
 
531   do 500 I=1,NN
        if(IC.eq.5) goto 549
        if(IC.eq.3)  Z=Z_min+(I-1)*Z_step
        if((IC.eq.1).or.(IC.eq.2).or.(IC.eq.4)) Xr=r_min+(I-1)*r_step
549       do 550 L=1,MS
             if(IC.eq.2) Z=Z_min+(L-1)*Z_step
             if((IC.eq.1).or.(IC.eq.3)) PHI=(PHI_min+(L-1)*
     1        PHI_step)*P/180
             if(IC.eq.5) then
                PHI=PHI5*P/180.
                Z=Z5
                Xr=r5
                goto 565
            endif
551       if(IC.eq.4) PHI=PHI4*P/180.
565       PHID=PHI*180./P
568       X(1)=Xr*COS(PHI)
          Y(1)=Xr*SIN(PHI)
cc
c     This part translates the coordinates ( PHIM ) to each magnetized piece
c  and calculates the advance of the easy axis on the rest of the pieces.
c
cc
          if (M.eq.1) goto 400
          do K=2,M
             PHIM=2.*P*(K-1)/M
             X(K)=X(1)*COS(PHIM)+Y(1)*SIN(PHIM)
             Y(K)=-X(1)*SIN(PHIM)+Y(1)*COS(PHIM)
          end do
cc
c  This part checks the computer's singularity points and adjusts .000001 in
c  the r dimension and .001 in ten thousands in the Z dimension to avoid
c  the chaos produced by the " smart man " --- computer or to avoid
c  the unrealistic value at the singularity points.
cc
400       do K=1,M
             if(X(K).eq.R) goto 433
             if(X(K).eq.(R+A)) goto 433
             if((X(K)-R)*SI34.eq.(Y(K)-B/2.)*CO34) goto 433
             if((X(K)-R-A)*SI34.eq.(Y(K)-D/2.)*CO34) goto 433
             if((X(K)-R)*SI34.eq.-(Y(K)+B/2.)*CO34) goto 433
             if((X(K)-R-A)*SI34.eq.-(Y(K)+D/2)*CO34) goto 433
           if(E.ne.0.) then
             if((X(K)-R-A)*SI56.eq.(Y(K)-D/2.)*CO56) goto 433
             if((X(K)-R-A-E)*SI56.eq.(Y(K)-F/2.)*CO56) goto 433
             if((X(K)-R-A)*SI56.eq.-(Y(K)+D/2.)*CO56) goto 433
             if((X(K)-R-A-E)*SI56.eq.-(Y(K)+F/2)*CO56) goto 433
           endif
          end do
          goto 436
433     if(Xr.lt.0) then
            Xr=Xr+.000001
           else
            Xr=Xr-.000001
          endif
          goto 568
436       if(Z.eq.C/2.) Z=Z+.001
          if(Z.eq.-C/2.) Z=Z-.001
cc
c     easy_axis is the angle that the easy axis advances.
cc    PHIM1 is the angle that a piece advances relative to the 1st piece.
          do II=1,M
             easy_axis=N*P*(II-1)/(M*1.)+AN
             PHIM1=2.*P/M*(II-1)
cc
c  This part calculates the atan functions.
cc
             OP1=OX1(X(II),Y(II),Z)
             OP2=OX2(X(II),Y(II),Z)
             OP3=OX3(X(II),Y(II),Z)
             OP4=OX4(X(II),Y(II),Z)
             OP5=OX5(X(II),Y(II),Z)
             OP6=OX6(X(II),Y(II),Z)
             OP7=OX7(X(II),Y(II),Z)
             OP8=OX8(X(II),Y(II),Z)
             OP9=OX9(X(II),Y(II),Z)
             OP10=OX10(X(II),Y(II),Z)
             OP11=OX11(X(II),Y(II),Z)
             OP12=OX12(X(II),Y(II),Z)
             OP13=OX13(X(II),Y(II),Z)
             OP14=OX14(X(II),Y(II),Z)
             OP15=OX15(X(II),Y(II),Z)
             OP16=OX16(X(II),Y(II),Z)
             TP1=OY1(X(II),Y(II),Z)
             TP2=OY2(X(II),Y(II),Z)
             TP3=OY3(X(II),Y(II),Z)
             TP4=OY4(X(II),Y(II),Z)
             TP5=OY5(X(II),Y(II),Z)
             TP6=OY6(X(II),Y(II),Z)
             TP7=OY7(X(II),Y(II),Z)
             TP8=OY8(X(II),Y(II),Z)
 
cc
c  This part calculates the log part
cc
             BYY1=BY1(X(II),Y(II),Z)
             BYY2=BY2(X(II),Y(II),Z)
             BYY3=BY3(X(II),Y(II),Z)
             BYY4=BY4(X(II),Y(II),Z)
             BYY5=BY5(X(II),Y(II),Z)
             BYY6=BY6(X(II),Y(II),Z)
             BYY7=BY7(X(II),Y(II),Z)
             BYY8=BY8(X(II),Y(II),Z)
             BZZ1=BZ1(X(II),Y(II),Z)
             BZZ2=BZ2(X(II),Y(II),Z)
             BZZ3=BZ3(X(II),Y(II),Z)
             BZZ4=BZ4(X(II),Y(II),Z)
             BZZ5=BZ5(X(II),Y(II),Z)
             BZZ6=BZ6(X(II),Y(II),Z)
             BZZ7=BZ7(X(II),Y(II),Z)
             BZZ8=BZ8(X(II),Y(II),Z)
             BXX1=BX1(X(II),Y(II),Z)
             BXX2=BX2(X(II),Y(II),Z)
             BXX3=BX3(X(II),Y(II),Z)
             BXX4=BX4(X(II),Y(II),Z)
             BXX5=BX5(X(II),Y(II),Z)
             BXX6=BX6(X(II),Y(II),Z)
             BXX7=BX7(X(II),Y(II),Z)
             BXX8=BX8(X(II),Y(II),Z)
             BL1=BLZ1(X(II),Y(II),Z)
             BL2=BLZ2(X(II),Y(II),Z)
             BL3=BLZ3(X(II),Y(II),Z)
             BL4=BLZ4(X(II),Y(II),Z)
             BL5=BLZ5(X(II),Y(II),Z)
             BL6=BLZ6(X(II),Y(II),Z)
             BL7=BLZ7(X(II),Y(II),Z)
             BL8=BLZ8(X(II),Y(II),Z)
             BTY=LOG(BYY7*BYY8*BYY3*BYY4/(BYY5*BYY6*BYY2*BYY1))
             BTZ=LOG(BZZ7*BZZ8*BZZ3*BZZ4/(BZZ5*BZZ6*BZZ2*BZZ1))
             BTX3=LOG(BXX1*BXX2/(BXX3*BXX4))
             BTX4=LOG(BXX5*BXX6/(BXX7*BXX8))
 
             BTZ3=LOG(BL3*BL4/(BL1*BL2))
             BTZ4=LOG(BL7*BL8/(BL5*BL6))
cc
c      contributions from surfaces 5 and 6
c
         if(E.ne.0.) then
             OP17=OX17(X(II),Y(II),Z)
             OP18=OX18(X(II),Y(II),Z)
             OP19=OX19(X(II),Y(II),Z)
             OP20=OX20(X(II),Y(II),Z)
             OP21=OX21(X(II),Y(II),Z)
             OP22=OX22(X(II),Y(II),Z)
             OP23=OX23(X(II),Y(II),Z)
             OP24=OX24(X(II),Y(II),Z)
             TP9=OY9(X(II),Y(II),Z)
             TP10=OY10(X(II),Y(II),Z)
             TP11=OY11(X(II),Y(II),Z)
             TP12=OY12(X(II),Y(II),Z)
             TP13=OY13(X(II),Y(II),Z)
             TP14=OY14(X(II),Y(II),Z)
             TP15=OY15(X(II),Y(II),Z)
             TP16=OY16(X(II),Y(II),Z)
             BXX9=BX9(X(II),Y(II),Z)
             BXX10=BX10(X(II),Y(II),Z)
             BXX11=BX11(X(II),Y(II),Z)
             BXX12=BX12(X(II),Y(II),Z)
             BL9=BLZ9(X(II),Y(II),Z)
             BL10=BLZ10(X(II),Y(II),Z)
             BL11=BLZ11(X(II),Y(II),Z)
             BL12=BLZ12(X(II),Y(II),Z)
             BL13=BLZ13(X(II),Y(II),Z)
             BL14=BLZ14(X(II),Y(II),Z)
             BL15=BLZ15(X(II),Y(II),Z)
             BL16=BLZ16(X(II),Y(II),Z)
             BTX5=LOG(BXX4*BXX10/(BXX1*BXX9))
             BTX6=LOG(BXX8*BXX12/(BXX5*BXX11))
             BTZ5=LOG(BL9*BL10/(BL11*BL12))
             BTZ6=LOG(BL15*BL16/(BL14*BL13))
             FX5=-(OP17+OP18+OP19+OP20)
             FX6=-(OP21+OP22+OP23+OP24 )
             FY5=TP9+TP10+TP11+TP12
             FY6=TP13+TP14+TP15+TP16
          endif
cc
c      summations and translations
c
212          FX1=-OP1-OP2-OP3-OP4+OP5+OP6+OP7+OP8
             FX3=-(OP9+OP10+OP11+OP12)
             FX4=-(OP13+OP14+OP15+OP16 )
 
             FY3=TP1+TP2+TP3+TP4
             FY4=TP5+TP6+TP7+TP8
 
 
             HXB=FX1*COS(easy_axis)+(FX4*SI34+CO34*BTX4)*SIN(easy_axis
     1         -theta34)-(FX3*SI34+CO34*BTX3)*SIN(easy_axis+theta34)
     1         -(FX5*SI56+CO56*BTX5)*SIN(easy_axis
     1         +theta56)+(FX6*SI56+CO56*BTX6)*SIN(easy_axis-theta56)
             HBX=HXB*BR(II)/4./P
             HYB=COS(easy_axis)*BTY+(FY4*CO34+SI34*BTX4)*SIN(easy_axis
     1         -theta34)-(FY3*CO34-SI34*BTX3)*SIN(easy_axis+theta34)
     1         +(FY6*CO56+SI56*BTX6)*SIN(easy_axis
     1         -theta56)-(FY5*CO56-SI56*BTX5)*SIN(easy_axis+theta56)
             HBY=HYB*BR(II)/4./P
             HX(II)=HBX*COS(PHIM1)-HBY*SIN(PHIM1)
             HY(II)=HBX*SIN(PHIM1)+HBY*COS(PHIM1)
             HZB=BTZ*COS(easy_axis)+BTZ4*SIN(easy_axis-theta34)+BTZ3*
     1         SIN(easy_axis+theta34)+BTZ6*SIN(easy_axis-theta56)-BTZ5*
     1         SIN(easy_axis+theta56)
             HZ(II)=BR(II)*HZB/4./P
          end do
          HHX=0.
          HHY=0.
          HHZ=0.
          do J=1,M
             HHX=HHX+HX(J)
             HHY=HHY+HY(J)
             HHZ=HHZ+HZ(J)
          end do
          BRA=HHX*COS(PHI)+HHY*SIN(PHI)
          BPH=-HHX*SIN(PHI)+HHY*COS(PHI)
c       The next 13 lines were added by H. Pfutzner at TUNL to allow 
c       including an additional solenoidal field
          BPERP=((BRA*BRA)+(BPH*BPH))**.5
          NZ=INT(Z+(C/2)+Z0+1.0)
          ZZ=Z+(C/2)+Z0+1.0-NZ
          BRAL=(Xr/2)*(I1*(BZL1(NZ+1)-BZL1(NZ))+
     +                 I2*(BZL2(NZ+1)-BZL2(NZ))+
     +                 I3*(BZL3(NZ+1)-BZL3(NZ)))
          HHZL=I1*(BZL1(NZ)+ZZ*(BZL1(NZ+1)-BZL1(NZ)))+
     +         I2*(BZL2(NZ)+ZZ*(BZL2(NZ+1)-BZL2(NZ)))+
     +         I3*(BZL3(NZ)+ZZ*(BZL3(NZ+1)-BZL3(NZ)))
          BRAT=BRA+BRAL
          HHZT=HHZ+HHZL
cc
c      output
cc
          BT=(BRA*BRA+BPH*BPH+HHZ*HHZ)**.5
          BTLH=(BRAT*BRAT+BPH*BPH+HHZT*HHZT)**.5
c
c         BRA is B radial from hexapole magnet
c         BPH is B azimuthal from hexapole magnet
c         HHZ is B axial from hexapole magnet
c         BPERP is vector sum of BRA and BPH
c         BT is vector sum of BRA, BPH, and HHZ
c
c         BRAL is B radial due to gradient in solenoidal magnetic field
c         HHZL is B axial from solenoidal magnet
c         BTLH is vector sum of BRA, BPH, HHZ, BRAL, and HHZL
c
          if(IC.eq.1) write(10,451)Xr,PHID,BRA,BPH,HHZ,BPERP,BT
     +                ,BRAL,HHZL,BTLH
          if(IC.eq.2) write(10,451)Xr,Z,BRA,BPH,HHZ,BPERP,BT
     +                ,BRAL,HHZL,BTLH
          if(IC.eq.3) write(10,451)Z,PHID,BRA,BPH,HHZ,BPERP,BT
     +                ,BRAL,HHZL,BTLH
          if(IC.eq.4) write(10,452)Xr,BRA,BPH,HHZ,BPERP,BT
     +                ,BRAL,HHZL,BTLH
          if(IC.eq.5) then
           write(10,455)Xr,PHID,Z,BRA,BPH,HHZ,BPERP,BT,BRAL,HHZL,BTLH
           goto 595
          endif
451       format (1x,2F8.2,5X,8F7.3)
452       format (1x,F8.2,13X,8F7.3)
455       format (1x,2F8.2,F5.2,8F7.3)
550       continue
500   continue
595   if(IC.eq.5) then
        type *,'  r=',Xr,'  phi=',PHID,'  z=',Z
        type *,' '
        type *,'  Br=',BRA,'   Bphi=',BPH,'   Bz=',HHZ
        type *,' '
        type *,'  Bt=',BT
      end if
 
596   type *,'  '
      type *,' Thank you, one more run? (Y/N)'
      read (5,FMT='(A)')REPLY
      ISTAT = STR$UPCASE(REPLY,REPLY)                    ! Set REPLY UPPER CASE
      if(REPLY.eq.'Y') then
         goto 600
       else
         goto 700
      end if
600   write(6,660)IN,M
660   format('   N = ',I2,'   M = ',I2,' ,   any change ? (Y/N)')
      read (5,fmt='(A)') CHANGE
      ISTAT = STR$UPCASE(CHANGE,CHANGE)                ! Set CHANGE UPPER CASE
      if(CHANGE.eq.'Y') goto 1
601   write(6,661)R,A,B,D,E,FF,C,AN1
661   format('  R =',f8.3,' , A =',f8.3,' , B =',f8.3,' , D =',f8.3,
     1/,'  E =',f8.3,' , F =',f8.3,' , C =',f10.2,' , AN1 =',f5.2) 
      type *,'    Any change ? (Y/N)'
      read (5,fmt='(A)') CHANGE
      ISTAT = STR$UPCASE(CHANGE,CHANGE)                  ! Set CHANGE UPPER CASE
      if(CHANGE.eq.'Y') goto 2
602   do k=1,M
         residual(1:7)='Brm(  )'
         write(residual(5:6),620)k
         write(6,621)residual,Br(k)
620      format(I2)
621      format(5x,A7,' = ',f6.3)
      end do
      type *, '   Any change ? (Y/N)'
      read (5,fmt='(A)') CHANGE
      ISTAT = STR$UPCASE(CHANGE,CHANGE)                 ! Set CHANGE UPPER CASE
      if(CHANGE.eq.'Y') goto 3
         goto 4
 
700   stop
      end