C
C PURPOSE: CALCULATES A MAGNETIC FIELD ACCORDING TO FONG, FAN 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
C COMMON:
C
C COMMENTS: AUGUST 24, 2000, P. JUDGE
C
C NOTES:    USER-SUPPLIED ROUTINE
C 
	SUBROUTINE FFL(RAD,THET,PH,BMAG,T,P,EDENS,HDENS,
     *    TEMP,VEL,TURBV,IPRINT)
	INCLUDE 'PREC'
	INCLUDE 'PARAM'
	PARAMETER(MAXD=10,MAXLEG=200)
        INCLUDE 'CCONST'
	INCLUDE 'CGRID'
	INCLUDE 'CINPUT'
	INCLUDE 'CLU'
	INCLUDE 'CINTS'
	COMMON /CFFL/ ACONST,B0,B1, ALPHAX, ALPHAY,P0,T0,ALEG(0:MAXLEG),NLEG
	COMMON /WATM/BX,BY,BZ
        DIMENSION PM(0:MAXLEG,0:MAXLEG),PD(0:MAXLEG,0:MAXLEG)
        DIMENSION PM2(0:MAXLEG,0:MAXLEG),PD2(0:MAXLEG,0:MAXLEG)
	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,'FFL.DAT')
           READ(LDUMS,*) ACONST
	   READ(LDUMS,*) B0
	   READ(LDUMS,*) B1
	   READ(LDUMS,*) ALPHAY, ALPHAX
	   READ(LDUMS,*) P0
	   READ(LDUMS,*) T0
	   READ(LDUMS,*) NLEG
	   IF(NLEG .GT. MAXLEG) CALL STOP('NLEG > MAXLEG')
	   READ(LDUMS,*)(ALEG(J),J=0,NLEG)
	   WRITE(LOUT,4001) ACONST,B0,B1,ALPHAY,ALPHAX,P0,T0,NLEG
 4001	   FORMAT('1 CORONAL ROUTINE IS FFL'/
     *     'A= ',F11.6,' B0= ',F10.4,' B1= ',F10.4,' ALPHAY= ',F10.4,
     *     ' ALPHAX= ',F10.4/'P0= ',F10.4,' T0=',1P,E10.3,0P/
     *     'NLEG= ',I4/'ALEG:')
	   WRITE(LOUT,4002)(ALEG(J),J=0,NLEG)
 4002	   FORMAT(1P,8(E9.2,1X),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    - THIS IS THE FORMALISM HANDED TO PGJ BY B.C. LOW NOVEMBER 14, 2002
C
	COSTH=Z1/R
	sign=1.
	if(z1 .lt. 0.) sign=-1.
	SINTH=SQRT(1.-COSTH*COSTH)
	TH=ATAN2(SINTH,COSTH)
	TANTH=SINTH/COSTH
C
	FF = (1.-R*R/ACONST/ACONST)
	GG = SQRT(FF*FF + 4.*R*R/ACONST/ACONST*COSTH*COSTH)
	XI =  SQRT(0.5*(-FF + GG))
	ETA = SQRT(0.5*( FF + GG))*SIGN
C	AT=ATAN(1./XI)
	AT=PI/2.0 - ATAN(XI) 
	BR0=B1*4./(ACONST*R*(XI*XI+ETA*ETA+EPSILON))
	BR1=-3.*XI*(5.*XI*XI+3.)*AT + 15.*XI*XI +4. 
	BR2=-3.*(1.+XI*XI)*(1.+5.*XI*XI)*AT + 15.*XI*XI*XI + 13.*XI
	BR= 2.*B0*COSTH/R/R/R - BR0*ETA*((1.-ETA*ETA)*(5.*ETA*ETA-1)*BR1 
     *      + XI*(3.-5.*ETA*ETA)*BR2)
        BTH= B0*SINTH/R/R/R - BR0*SQRT((1.-ETA*ETA)/SQRT(1.+XI*XI))
     *       * (XI*SQRT(1.+XI*XI)*(5.*ETA*ETA-1.)*BR1
     *       - ETA*ETA*(3.-5.*ETA*ETA)*BR2)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  ADD IN LEGENDRE EXPANSION OF IMAGE POTENTIAL
C
	IONE=1
	CALL PNL(MAXLEG,IONE,NLEG,COSTH,PM,PD)
	SUMR=ZERO
	SUMTH=ZERO
	RVAR=ONE/R/R
	DO 101 N=1,NLEG
           RVAR=RVAR/R
	   DELR=ALEG(N)*RVAR*(PM(1,N)/TANTH-PD(1,N)*SINTH)
	   DELTH=ALEG(N)*N*RVAR*PM(1,N)*SINTH
	   SUMR=SUMR+DELR
	   SUMTH=SUMTH+DELTH
 101	CONTINUE
C	IF((0.99 .LT. R) .AND. (1.01 .GT. R)) THEN 
C     	   WRITE(LOUT,2007) BR, SUMR,DELR/SUMR
C 2007	   FORMAT('BR SUMR DELR/SUMR'/3(F11.5,1X))
C	ENDIF
	BR=BR+SUMR
 	BTH=BTH+SUMTH
C
C END OF LEGENDRE EXPANSION
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
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