C
C PURPOSE: CALCULATES THE S-FRAME IRRED. TENSOR COMPONENTS OF AVG RADN FIELD
C
C  THIS SUBROUTINE CALCULATES THE IRREDUCIBLE TENSOR COMPONENTS
C  OF THE AVERAGE RADIATION FIELD FROM THE PHOTOSPHERE IN THE
C  S FRAME, GIVEN THE LIMB-DARKENING FUNCTION, AND ASSUMING
C  CYLINDRICAL SYMMETRY OF THE RADIATION FIELD. 
C
C  AN 8-POINT GAUSSIAN QUADRATURE FORMULA IS ADOPTED FOR THE 
C  CALCULATION OF THE WEIGHTED INTEGRALS. THE INTEGRAL IN 
C  (Y;AMU_M,1) IS TRANSFORMED INTO AN INTEGRAL IN (X;-1,1) FOR 
C  THE PURPOSE OF GAUSSIAN INTEGRATION
C
C  H IS THE HEIGHT OF THE OBSERVED POINT OVER THE LIMB, EXPRESSED
C  IN SOLAR RADII
C
C INPUTS:
C       KR    INDEX OF THE TRANSITION
C       (FLIMB EXTERNAL  FUNCTION THAT HAS TO BE DEFINED)
C
C OUTPUTS:
C       RJ0   K=0 TENSOR COMPONENT OF THE RADIATION FIELD
C       RJ1   K=1 TENSOR COMPONENT OF THE RADIATION FIELD
C                 (THIS IS ZERO IN CYLINDRICAL SYMMETRY)
C       RJ2   K=2 TENSOR COMPONENT OF THE RADIATION FIELD
C             SEE EQUATIONS 22 AND 23 OF 
C             CASINI, R. & JUDGE, P. G., 1999. AP J 522, 524-539
C
C COMMON:
C       CORON
C COMMENTS: OCTOBER 6, 1999, P. JUDGE
C 
	SUBROUTINE FIELD_INT(KR,RJ0,RJ1,RJ2)
	INCLUDE 'PREC'
	INCLUDE 'CORON'
	INCLUDE 'CINTS'
	DIMENSION X(4),W(4)
	DATA X/.183434642495650,.525532409916329,
     /	 .796666477413627,.960289856497536/
	DATA W/.362683783378362,.313706645877887,
     /	 .222381034453374,.101228536290376/
	EXTERNAL FLIMB
C  STATEMENT FUNCTION
	T2_0(X1)=(.5/SQRT(TWO))*(THREE*X1*X1-ONE)
	AMU_M=SQRT(H*(H+TWO))/(H+ONE)
	FACTOR1=.5*(ONE-AMU_M)			
	FACTOR2=FACTOR1+AMU_M
C  INITIALIZE
	RJ0=ZERO
	RJ2=ZERO
	DO I=1,4
	  YPOS=FACTOR1*X(I)+FACTOR2
	  YNEG=-YPOS+TWO*FACTOR2
	  TMP0POS=FLIMB(KR,YPOS,AMU_M)
	  TMP2POS=T2_0(YPOS)*TMP0POS
	  TMP0NEG=FLIMB(KR,YNEG,AMU_M)
	  TMP2NEG=T2_0(YNEG)*TMP0NEG
	  RJ0=W(I)*(TMP0POS+TMP0NEG)+RJ0
	  RJ2=W(I)*(TMP2POS+TMP2NEG)+RJ2
	END DO
	RJ0=.5*FACTOR1*RJ0
	RJ1=ZERO
	RJ2=.5*FACTOR1*RJ2
	RETURN
	END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C
C PURPOSE: CALCULATES RADIATION FIELD INTENSITY IN THE DIRECTION THETA=ARCOS(X)
C
C INPUTS:
C       KR  TRANSITION INDEX
C	X   COS THETA
C       XM  COS THETAM
C           YMU IS MU= COS ZETA IN EQUATION 29 OF CASINI & JUDGE 1999
C OUTPUTS:
C       FLIMB
C COMMON:
C
C RESTRICTIONS:
C       
C COMMENTS: OCTOBER 6, 1999, P. JUDGE
C 
	FUNCTION FLIMB(KR,X,XM)
	INCLUDE 'PREC'
	INCLUDE 'PARAM'
	INCLUDE 'CATOM'
	INCLUDE 'CCONST'
	INCLUDE 'CINTS'
	YMU=SQRT((X+XM)*(X-XM)/(ONE-XM*XM))
C
C  LAMBDA IS IN ANGSTROMS
C  XNU IS ANGULAR FREQUENCY OF TRANSITION
C
	WL=ALAMB(KR)
	WM=WL/1.E4
	XNU=1.E8*(TWO*PI)*(CC/WL)
C
C XT IS DISK CENTER RADIATION TEMPERATURE FOR I (ALLEN 1973, SEC 80)
C
	XT=6050.
C
C ALLEN'S (1973) LIMB DARKENING (P 170, TOP)
C
	CALL UV(WM,U,V)
	FLIMB=PLANCK(XNU,XT)*(ONE-U-V+(U+V*YMU)*YMU)
	RETURN
	END