C
C PURPOSE: CALCULATES A MAGNETIC FIELD ACCORDING TO ZHANG AND LOW
C
C INPUTS:
C       RAD    RADIUS OF POINT RELATIVE TO SUN CENTER, UNITS OF SOLAR RADII
C       THET   AZIMUTHAL ANGLE, UNITS RADIANS, IN S FRAME
C       PH     POLAR ANGLE, UNITS RADIANS, IN S FRAME
C       IPRINT SET TO 1 FOR MORE VERBOSE OUTPUT
C OUTPUTS:
C       B      MAGNITUDE OF FIELD
C       T      POLAR ANGLE
C       P      AZIMUTHAL ANGLE
C       EDENS  ELECTRON DENSITY /CM3
C       HDENS  HYDROGEN NUCLEI DENSITY /CM3 
C       TEMP   ELECTRON TEMPERATURE K
C       VEL    FLOW VELOCITY ALONG LINE OF SIGHT (+VE MEANS TOWARDS OBSERVER)
C       TURBV  MICROTURBULENCE KM/S (CONTRIBUTES TO LINE PROFILE WIDTHS)
C COMMON:
C
C COMMENTS: APRIL 19, 2006
C
C NOTES:    USER-SUPPLIED ROUTINE
C 
	SUBROUTINE ZL(RAD,THET,PH,BMAG,T,P,EDENS,HDENS,
     *    TEMP,VEL,TURBV,IPRINT)
	INCLUDE 'PREC'
	INCLUDE 'PARAM'
        INCLUDE 'CCONST'
	INCLUDE 'CGRID'
	INCLUDE 'CINPUT'
	INCLUDE 'CLU'
	INCLUDE 'CINTS'
	COMMON /CZL/ ACONST,B0,CI0, ALPHAX, ALPHAY,P0,T0
	COMMON /WATM/BX,BY,BZ
	DATA KCALL /0/
	DATA EPSILON, ZER, SMALL /1.E-3,0.,1.E-20/
	SAVE KCALL
C
	IF (KCALL .EQ. 0) THEN 
	   CALL CSTRIP(LDIP,'ZL.DAT')
           READ(LDUMS,*) ACONST
	   READ(LDUMS,*) B0
	   READ(LDUMS,*) CI0
	   READ(LDUMS,*) ALPHAY, ALPHAX
	   READ(LDUMS,*) P0
	   READ(LDUMS,*) T0
	   WRITE(LOUT,4001) ACONST,B0,CI0,ALPHAY,ALPHAX,P0,T0
 4001	   FORMAT('1 CORONAL ROUTINE IS ZL'/
     *     'A= ',F11.6,' B0= ',F10.4,' CI0= ',F10.4,' ALPHAY= ',F10.4,
     *     ' ALPHAX= ',F10.4/'P0= ',F10.4,' T0=',1P,E10.3,0P)
	   CALL CLOSE(LDUMS)
	   KCALL=KCALL+1
	ENDIF
C
C  CALCULATE CARTESIAN COORDINATES OF THE POINT RAD,THET,PH
C  X IS LOS
C
	X=RAD*SIN(THET)*COS(PH)
	Y=RAD*SIN(THET)*SIN(PH)
	Z=RAD*COS(THET)
        R=RAD
C
C ROTATE TO (R,THETA) PLANE IN WHICH FFL SOLUTION IS GIVEN IN TERMS
C OF STREAM FUNCTION. 
C	
C ROTATION ABOUT Y AXIS, AND THEN X-AXIS
C
	X1=X
	Y1=Y
	Z1=Z
	CALL ROT(X1,Z1,ALPHAY)
	CALL ROT(Y1,Z1,ALPHAX)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C ROTATED FRAME:
C  CALCULATE THE (R,THETA) COMPONENTS OF THE MAGNETIC FIELD IN THE ROTATED 
C  REFERENCE FRAME. 
C
	COSTH=Z1/R
	SINTH=SQRT(1.-COSTH*COSTH)
	TH=ATAN2(SINTH,COSTH)
	TANTH=SINTH/COSTH
C
	DTH=0.01
	DR=0.001*R
C
C SIMPLE CENTERED DIFFERENCE, 4 CALLS
C
C	BR=1./R/R/SINTH*(AZL(R,TH+DTH,CI0,ACONST,B0) - 
C     *     AZL(R,TH-DTH,CI0,ACONST,B0))/2./DTH
C	BTH=-1./R/SINTH * (AZL(R+DR,TH,CI0,ACONST,B0) - 
C     *     AZL(R-DR,TH,CI0,ACONST,B0))/2./DR
C ALL DONE!
C
C SIMPLE CENTERED DIFFERENCE, NR ITERATION
C
	BR=1./R/R/SINTH*DAZLDTH(R,TH,CI0,ACONST,B0)
	BTH=-1./R/SINTH*DAZLDR(R,TH,CI0,ACONST,B0)
C 
C CARTESIAN COMPONENTS IN THIS ROTATED REFERENCE FRAME
C
	BMAG=SQRT(BR*BR+BTH*BTH)
	V=SQRT(X1*X1+Y1*Y1)
	COSPH=X1/V
	SINPH=Y1/V
        BZ=  BR*COSTH - BTH*SINTH
	BX= (BR*SINTH + BTH*COSTH)*COSPH
	BY= (BR*SINTH + BTH*COSTH)*SINPH
C
C DENSITY RESTRICTION AROUND CURRENT SHEET, NARROWING WITH R
C NOT IN PRESSURE EQUILIBRIUM 
C BEGIN
C	THICK=0.1/R/R
C	TH1=ACOS(COSTH)
C       DELTA=(PI/2.-TH1)/THICK
C END
C ROTATED FRAME END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C TRANSFORM BACK, ROTATE ABOUT X AND Y AXES
C
	BX1=BX
	BY1=BY
	BZ1=BZ
	CALL ROT(BY1,BZ1,-ALPHAX)
	CALL ROT(BX1,BZ1,-ALPHAY)
C  RE-ASSIGN
	BX=BX1
	BY=BY1
	BZ=BZ1
C
C  RETURN VALUES OF B,T,P
C
	BMAG=SQRT(BX*BX+BY*BY+BZ*BZ)
	T=ACOS(BZ/BMAG)
	P=ATAN2(BY,BX)
C
C THERMAL PARAMETERS
C
C
C  ELECTRON DENSITY AND TEMPERATURE.
C
	VEL=0.
	TURBV=15.
C	DELTA=0.
C	GFAC=EXP(-DELTA*DELTA)
	GRAV=2.7398E4
	TEMP=T0
C	EDENS=P0/BK/T0/1.88*GFAC
	EDENS=P0/BK/T0/1.88
	RHO=1.3625*0.8*UU
	PTOT=BK*TEMP*1.88 + 0.5* RHO*(TURBV*TURBV)*1.E10
	HSCALE = PTOT/RHO/GRAV
	HSCALE=HSCALE/RSUNCM
	FFAC= MAX((R-1.),ZER)
	FFAC=EXP(-FFAC/HSCALE)
	EDENS=EDENS*FFAC
	HDENS=0.8*EDENS
	RETURN
	END
C
C**********************************************************************
C
	FUNCTION AZL(R, THETA, CI0, A0, B0)
	INCLUDE 'PREC'
	INCLUDE 'CCONST'
	EPS = 1.E-6
	AZL=0.0
	ST=SIN(THETA)
C   
C ABRAMOWITZ & STEGUN  EQS. 25.4.38, INTEGRAL OF F(X)/SQRT(1-X^2) DX BETWEEN
C -1 AND 1
C HERE, F(X) IS X/SQRT(A^2+R^2 - 2AR SIN(TH) SQRT(1-X^2))
C
	NMAX=40
	SUM=0.
	W = PI/NMAX
	DO 100 I=1,NMAX
	   X = COS( PI*(2.*I-1.)/NMAX )
	   F = X/SQRT(EPS+A0*A0+R*R -2.*A0*R*ST*X)
	   FIMAGE = X/SQRT(EPS+A0*A0*R*R +1. -2.*A0*R*ST*X)
	   SUM=SUM+F-FIMAGE
 100	CONTINUE
	AZL = 2.*CI0*A0*R*ST*SUM*W + B0*ST*ST/R
	RETURN
	END
C
C**********************************************************************
C BASED ON DFRIDR.F OF NUMERICAL RECIPES
C
      FUNCTION DAZLDTH(R,TH,CI0,ACONST,B0)
      INTEGER NTAB
      REAL DAZLDTH,ERR,H,X,FUNC,CON,CON2,BIG,SAFE
      PARAMETER (CON=1.4,CON2=CON*CON,BIG=1.E30,NTAB=10,SAFE=2.)
      INTEGER I,J
      REAL ERRT,FAC,HH,A(NTAB,NTAB)
      X=TH
      H=0.01
      HH=H
      A(1,1)=(AZL(R,TH+HH,CI0,ACONST,B0)-AZL(R,TH-HH,CI0,ACONST,B0))
     * /(2.0*HH)
      ERR=BIG
      DO 12 I=2,NTAB
        HH=HH/CON
        A(1,I)=(AZL(R,TH+HH,CI0,ACONST,B0)-AZL(R,TH-HH,CI0,ACONST,B0))
     *   /(2.0*HH)
        FAC=CON2
        DO 11 J=2,I
          A(J,I)=(A(J-1,I)*FAC-A(J-1,I-1))/(FAC-1.)
          FAC=CON2*FAC
          ERRT=MAX(ABS(A(J,I)-A(J-1,I)),ABS(A(J,I)-A(J-1,I-1)))
          IF (ERRT.LE.ERR) THEN
            ERR=ERRT
            DAZLDTH=A(J,I)
          ENDIF
11      CONTINUE
        IF(ABS(A(I,I)-A(I-1,I-1)).GE.SAFE*ERR) RETURN
12    CONTINUE
      RETURN
      END
C
C**********************************************************************
C
      FUNCTION DAZLDR(R,TH,CI0,ACONST,B0)
      INTEGER NTAB
      REAL DAZLDTH,ERR,H,X,FUNC,CON,CON2,BIG,SAFE
      PARAMETER (CON=1.4,CON2=CON*CON,BIG=1.E30,NTAB=10,SAFE=2.)
      INTEGER I,J
      REAL ERRT,FAC,HH,A(NTAB,NTAB)
      X=R
      H=0.01*R
      HH=H
      A(1,1)=(AZL(R+HH,TH,CI0,ACONST,B0)-AZL(R-HH,TH,CI0,ACONST,B0))
     * /(2.0*HH)
      ERR=BIG
      DO 12 I=2,NTAB
        HH=HH/CON
        A(1,I)=(AZL(R+HH,TH,CI0,ACONST,B0)-AZL(R-HH,TH,CI0,ACONST,B0))
     *   /(2.0*HH)
        FAC=CON2
        DO 11 J=2,I
          A(J,I)=(A(J-1,I)*FAC-A(J-1,I-1))/(FAC-1.)
          FAC=CON2*FAC
          ERRT=MAX(ABS(A(J,I)-A(J-1,I)),ABS(A(J,I)-A(J-1,I-1)))
          IF (ERRT.LE.ERR) THEN
            ERR=ERRT
            DAZLDR=A(J,I)
          ENDIF
11      CONTINUE
        IF(ABS(A(I,I)-A(I-1,I-1)).GE.SAFE*ERR)RETURN
12    CONTINUE
      RETURN
      END
