C
C PURPOSE: CALCULATES A SIMPLE DIPOLAR FIELD
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: AUGUST 24, 2000, P. JUDGE
C
C NOTES:    USER-SUPPLIED ROUTINE
C 
	SUBROUTINE DIPOLE(RAD,THET,PH,B,T,P,EDENS,HDENS,
     *    TEMP,VEL,TURBV,IPRINT)
	INCLUDE 'PREC'
	PARAMETER(MAXD=10)
        INCLUDE 'CCONST'
	INCLUDE 'CGRID'
	INCLUDE 'CINPUT'
	INCLUDE 'CLU'
	COMMON /CDIP/ XD,YD,ZD,XM,YM,ZM,DM,DELTA,BASE,NTUBE
        COMMON /ADIP/ AALP(MAXD),DALP(MAXD),AP(MAXD),DP(MAXD)
	COMMON /WATM/BX,BY,BZ
	DATA ICALL, ZER, SMALL /0,0.,1.E-20/
	SAVE ICALL
C
	IF (ICALL .EQ. 0) THEN 
	   CALL CSTRIP(LDIP,'DIPOLE.DAT')
	   READ(LDUMS,*) RD,THD,PHID
	   THD=THD*PI/180.
	   PHID=PHID*PI/180.
	   XD=RD*SIN(THD)*COS(PHID)
	   YD=RD*SIN(THD)*SIN(PHID)
	   ZD=RD*COS(THD)
	   READ(LDUMS,*) XM,YM,ZM
	   READ(LDUMS,*) DM
	   READ(LDUMS,*) DELTA
	   READ(LDUMS,*) BASE
	   READ(LDUMS,*) NTUBE
	   IF(NTUBE .GT. MAXD) CALL STOP('NTUBE .GT. MAXD') 
	   DO I=1,NTUBE
	      READ(LDUMS,*) DAA,DBB,DCC,DDD
	      AALP(I)=DAA
	      DALP(I)=DBB
	      AP(I)=DCC
	      DP(I)=DDD
	   END DO
	   CALL CLOSE(LDUMS)
	   ICALL=ICALL+1
	ENDIF
C
C   NORMALIZE TO MODULO 1
C
	TOT=SQRT(XM*XM+YM*YM+ZM*ZM)
	XM=XM/TOT
	YM=YM/TOT
	ZM=ZM/TOT
C
C  CALCULATE CARTESIAN COORDINATES OF THE POINT RAD,THET,PH
C
	X=RAD*SIN(THET)*COS(PH)
	Y=RAD*SIN(THET)*SIN(PH)
	Z=RAD*COS(THET)
        RR=SQRT(X*X+Y*Y+Z*Z)
C
C  CALCULATE THE SCALAR MAGNETIC POTENTIAL
C  DM IS DIPOLE MOMENT
C
	RX=X-XD
	RY=Y-YD
	RZ=Z-ZD
C	V=DIPOT(DM,XM,YM,ZM,RX,RY,RZ)
C
C  CALCULATE MAGNETIC FIELD AT RAD, THET, PH
C  NUMERICAL DERIVATIVES ARE USED (CENTERED)
C
	DX=DELTA
	VP1=DIPOT(DM,XM,YM,ZM,RX + DX,RY,RZ) 
	VM1=DIPOT(DM,XM,YM,ZM,RX - DX,RY,RZ) 
	BX= 0.5*(VP1-VM1)/DX
C
	DY=DELTA
	VP1=DIPOT(DM,XM,YM,ZM,RX,RY+DY,RZ) 
	VM1=DIPOT(DM,XM,YM,ZM,RX,RY-DY,RZ) 
	BY= 0.5*(VP1-VM1)/DY
C
	DZ=DELTA
	VP1=DIPOT(DM,XM,YM,ZM,RX,RY,RZ+DZ) 
	VM1=DIPOT(DM,XM,YM,ZM,RX,RY,RZ-DZ) 
	BZ= 0.5*(VP1-VM1)/DZ
C
	B=SQRT(BX*BX+BY*BY+BZ*BZ)
	T=ACOS(BZ/B)
	P=ATAN2(BY,BX)
C
C  ELECTRON DENSITY AND TEMPERATURE.
C  FILL UP A "CONSTANT ALPHA" FLUX TUBE WITH PLASMA
C  THE FLUX TUBE IS FLAGGED IN TWO WAYS:
C   1.  A VALUE OF ALPHA IS CHOSEN (THIS SPECIFIES THE 
C       "ONION SKIN" IN THE CYLINDRICALLY SYMMETRIC FIELD, AND 
C   2.  A VALUE FOR THE AZIMUTHAL ANGLE IS CHOSEN. THE DENSITY IS THEN
C       LOADED ON TO THE FLUX TUBE WITH A GAUSSIAN DISTRIBUTION CENTERED
C       ON THESE VALUES.
C
C  PP IS THE AZIMUTHAL COORDINATE IN THE REF FRAME OF THE DIPOLE
C  MUST FIRST CALCULATE COORDINATES OF THE POINT IN THE DIPOLAR REST FRAME
C
	BET=ATAN2(YM,ZM)
	ALP=-ATAN2(XM,SQRT(ZM*ZM+YM*YM))	
	RRX=COS(ALP)*RX + SIN(ALP)*SIN(BET)*RY + SIN(ALP)*COS(BET)*RZ
	RRY=COS(BET)*RY-SIN(BET)*RZ
	RRZ=-SIN(ALP)*RX + COS(ALP)*SIN(BET)*RY +COS(ALP)*COS(BET)*RZ
	PP=ATAN2(RRY,RRX)
C
C ZZ IS THE "Z COORDINATE" IN THE REF FRAME OF THE DIPOLE. 
C
	ZZ=RRZ
	R2=RRX*RRX + RRY*RRY + RRZ*RRZ
	ALPHA = DM/SQRT(R2)*(1.-ZZ*ZZ/R2)
C
C GAUSSIANS
C
	DO I=1,NTUBE
	   DA0=(ALPHA-AALP(I))/DALP(I)
	   DP0=(PP-AP(I))/DP(I)
	   EXPON=-DA0*DA0-DP0*DP0
	   IF(I .EQ. 1) THEN 
	      EDENS=MAX(BASE*EXP(EXPON),SMALL)
	   ELSE
	      EDENS=MAX(EDENS,BASE*EXP(EXPON))
	   ENDIF
	END DO
C
C HARD EDGES
C
C        EDENS=1.
C        IF(ABS(DA0) .LT. 1. .AND. ABS(DP0) .LT. 1.) EDENS=BASE
	TEMP=2.E6
	VEL=0.
	TURBV=20.
	GRAV=2.7398E4
	TEMP1=2.E6
	RHO=(1.3625*0.8*UU)
	PTOT=BK*TEMP1*1.88 + 0.5* RHO*(TURBV*TURBV)*1.E10
	HSCALE = PTOT/RHO /GRAV
	HSCALE=HSCALE/6.9599E10
	FFAC= MAX((RR-1.),ZER)
	FFAC=EXP(-FFAC/HSCALE)
C	EDENS=EDENS*FFAC
	HDENS=0.8*EDENS
C	IF(IWATMO .NE. 0 .AND. IPRINT .NE. 0) THEN 
C	   WRITE(*,*)' '
C	   WRITE(*,200)'DIPOLE: RAD THET PH ', RAD, THET, PH
C	   WRITE(*,200)'DIPOLE:  X  Y  Z      ',  X,  Y,  Z
C	   WRITE(*,200)'DIPOLE: RX RY RZ      ', RX, RY, RZ
C	   WRITE(*,200)'DIPOLE: XM YM ZM      ', XM, YM, ZM
C	   WRITE(*,200)'DIPOLE: RRX RRY RRZ   ', RRX, RRY, RRZ
C	   WRITE(*,200)'DIPOLE: VM1 VP1 V     ',VM1,VP1, V
C	   WRITE(*,200)'DIPOLE: BX BY BZ      ', BX, BY, BZ
C	   WRITE(*,200)'DIPOLE: B T P         ', B, T, P
C	   WRITE(*,201)'DIPOLE: EDENS,ALPHA,PP,DA,DP,(RR-1)/H  ', 
C     *       LOG10(EDENS),ALPHA,PP,DA0,DP0,(RR-1.)/HSCALE
C	   WRITE(*,200)'DIPOLE: ALP,BET,PP    ', ALP,BET,PP
C	ENDIF
C 200	FORMAT(A,3(2X,F8.3,1X))
C 201	FORMAT(A,6(1X,F7.2))
	RETURN
	END
C
C**********************************************************************
C
	FUNCTION DIPOT(DM,XM,YM,ZM,RX,RY,RZ)
	INCLUDE 'PREC'
C
	R=SQRT(RX*RX + RY*RY + RZ*RZ)
	DIPOT = DM*(XM*RX+YM*RY+ZM*RZ)/R/R/R
	RETURN
	END