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 OUTPUTS:
C       B      MAGNITUDE OF FIELD
C       T      POLAR ANGLE
C       P      AZIMUTHAL ANGLE
C              - THERMAL PARAMETERS
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)
        INCLUDE 'CCONST'
	INCLUDE 'CGRID'
	INCLUDE 'CINPUT'
	INCLUDE 'CLU'
	COMMON /CFFL/ B0,B1, ALPHAX, ALPHAY,P0,T0
	COMMON /WATM/BX,BY,BZ
	DATA KCALL /0/
	DATA EPSILON, ACONST, ZER, SMALL /1.E-3,1.1180340,0.,1.E-20/
	SAVE KCALL
C
	IF (KCALL .EQ. 0) THEN 
	   CALL CSTRIP(LDIP,'FFL.DAT')
	   READ(LDUMS,*) B0
	   READ(LDUMS,*) B1
	   READ(LDUMS,*) ALPHAY, ALPHAX
	   READ(LDUMS,*) P0
	   READ(LDUMS,*) T0
	   CALL CLOSE(LDUMS)
	   KCALL=KCALL+1
	   WRITE(*,*) 'CHECK FFL', 'B0,B1,ALPHAY,ALPHAX,P0,T0'
	   WRITE(*,*) 'CHECK FFL', B0,B1,ALPHAY,ALPHAX,P0,T0
	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
	SINTH=SQRT(1.-COSTH*COSTH)
	TH=ATAN2(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))
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)
	BMAG=SQRT(BR*BR+BTH*BTH)
C 
C CARTESIAN COMPONENTS IN THIS ROTATED REFERENCE FRAME
C
	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
 1003	format(1p,3(e9.2,1x),A,0p)
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