C
C PURPOSE: CALCULATE EMISSION COEFFICIENTS FOR STOKES VECTOR 
C
C INPUTS:
C       KR   TRANSITION INDEX (1=FIRST TRANSITION, ETC)
C       T0   GEOMETRIC TENSOR (TABLE 1 OF
C            CASINI, R. & JUDGE, P. G., 1999. AP J 522, 524-539)
C OUTPUTS:
C   EMISS(KR,I,NY)  EMISSION COEFFICIENT IN PH/CM2/S/SR/ANGSTROM 
C   WHERE  KR IS THE LINE INDEX
C          I IS THE STOKES INDEX (0=INTENSITY, 1=Q, ETC.)
C          NU IS THE FREQUENCY INDEX
C
C COMMON:
C       CSE CSLINE CATOM CORON SGNM CLU
C COMMENTS: OCTOBER 6, 1999, P. JUDGE
C 
C       4TH NOVEMBER 2019
C	ADDED WHITE LIGHT EMISSION TO EMISSIVITY
C	22 JULY 2022 COMMENTED OUT UNTIL FULLY TESTED
C
C	JULY 1 2022
C	PGJ FIX PROFILE NORMALIZATION BUG
C	THE ORIGINAL VERSION RETURNED IQUV VALUES
C       LARGER BY AN ERRONEOUS FACTOR OF 3.155 * (QNORM/10)	
C
	SUBROUTINE EMISSION(KR,T0)
	INCLUDE 'PREC'
	INCLUDE 'PARAM'
	INCLUDE 'CSE'
	INCLUDE 'CSLINE'
	INCLUDE 'CCONST'
	INCLUDE 'CATOM'
	INCLUDE 'CORON'
	INCLUDE 'SGNM'
	INCLUDE 'CINTS'
	INCLUDE 'CLU'
C
	PARAMETER (SQRT2=.14142135623730950488E1)
	PARAMETER (SQRT3=.17320508075688772935E1)
C
	DIMENSION T0(0:3,0:2)
	DIMENSION EM_FACT(0:4)
C
	IJ0=IRAD(KR)
	IJ1=JRAD(KR)
	RJ0=AJ(IJ0)
	RJ1=AJ(IJ1)
C
C CALCULATE C_COEFF
C
	HNU=HH*CC / (ALAMB(KR)/1.E8)
	C_COEFF=TOTN*SQRT(TWO*RJ1+ONE)*RHO(IJ1,0)*A(KR)
     *         *HNU/FOUR/PI
C
C       diagnostic output
C
C 	WRITE(*,*) 'KR,IJ1,C_COEFF,TOTN,RJ1,RHO(IJ1,0),A(KR)',
C     *     KR,IJ1,C_COEFF,TOTN,RJ1,RHO(IJ1,0),A(KR)
C
C  	write(*,*) kr, alamb(kr),label(ij1),rho(ij1,0)
	G0=GLAND(IJ0)
	G1=GLAND(IJ1)
	GEFF=0.5*(G0+G1)+0.25*(G0-G1)*(RJ0*(RJ0+ONE)-RJ1*(RJ1+ONE))
	SIGMA=RHO(IJ1,2)/RHO(IJ1,0)
C
C       diagnostic output
C        WRITE(*,*) RHO(IJ1,2),RHO(IJ1,0)
C
C   CALCULATE DCOEFF AND ECOEFF
C
	D_COEFF=-SG(NINT(RJ0+RJ1))*SQRT(THREE*(TWO*RJ1+1))
     /	 *FUN6J(ONE,ONE,TWO,RJ1,RJ1,RJ0)
C
C ORIGINAL 2006 CODE
C      E_COEFF=THREE*SQRT(TWO*RJ1+ONE)
C     /	 *(SG(NINT(RJ1-RJ0))*G1*SQRT(RJ1*(RJ1+ONE)*(TWO*RJ1+ONE))
C    /	 *FUN6J(ONE,TWO,ONE,RJ1,RJ1,RJ1)
C     /	 *FUN6J(ONE,ONE,ONE,RJ1,RJ1,RJ0)
C     /	 -G0*SQRT(RJ0*(RJ0+ONE)*(TWO*RJ0+ONE))
C     /	 *FUN9J(ONE,TWO,ONE,RJ0,RJ1,ONE,RJ0,RJ1,ONE))
c
c       ALTERNATE CALCULATION OF E GIVEN D AND GEFF
c	
	F_COEFF=3./4. * (RJ1*(RJ1+ONE) - RJ0*(RJ0+ONE) -2.)*(G1-G0)
	E_COEFF=D_COEFF/SQRT2*(F_COEFF+GEFF)
C
C   CALCULATE EMISSION COEFFICIENTS
C
C       PGJ 9 JUNE 2022 BEGIN
C      OLD;
C      ALARMOR=(1.0E-8*(0.25E0/PI)*(EE/EM)*BFIELD)/QNORM
C
C       NEW:
C       FLARMOR LARMOR FREQUENCY UNITS OF HZ
C
	FLARMOR = 1.E-8 * EE * BFIELD / EM / 4./ PI
C
C       PGJ 9 JUNE 2022 END
C	
C       LINE WIDTH IN DOPPLER UNITS CM/S
C	
	WCGS=DNYD*QNORM*1.E5
C
C       IN FREQUENCY UNITS HZ
C
	WHZ = WCGS * 1.E8 / ALAMB(KR)
	EM_FACT(0)=ONE+D_COEFF*SIGMA*T0(0,2)
	EM_FACT(1)=D_COEFF*SIGMA*T0(1,2)
	EM_FACT(2)=D_COEFF*SIGMA*T0(2,2)
	EM_FACT(3)=-(SQRT2/SQRT3)*FLARMOR/WHZ * (GEFF+E_COEFF*SIGMA)*T0(3,1) 
	EM_FACT(4)=-(SQRT2/SQRT3)*FLARMOR/WHZ *  GEFF*T0(3,1)
C       diagnostic output
C	write(*,*) em_fact(0),em_fact(1),em_fact(2),em_fact(3),em_fact(4)
C
C EMISS(*,4,*): MAGNETOGRAM STOKES V
C
C	SI=0.
C	SQ=0.
C	SU=0.
C	
C       4TH NOVEMBER 2019
C	CALL KCOR(KR, SI, SQ, SU)
C       4TH NOVEMBER 2019
C  	SG: note-- SCL below not in PGJ version
C         I am keeping in (commented)
C          because it is in commented lines below
C         Although no longer used in EMISS definitions
C	SCL=C_COEFF*1.E-13*ALAMB(KR)
C	DO NY=1,NQ(KR)
C	  EMISS(KR,0,NY)=EM_FACT(0)*PHI(KR,NY)*C_COEFF + SI
C	  EMISS(KR,1,NY)=EM_FACT(1)*PHI(KR,NY)*C_COEFF + SQ
C	  EMISS(KR,2,NY)=EM_FACT(2)*PHI(KR,NY)*C_COEFF + SU
C	  EMISS(KR,3,NY)=-EM_FACT(3)*PHIP(KR,NY)*SCL
C	  EMISS(KR,4,NY)=-EM_FACT(4)*PHIP(KR,NY)*SCL
C	END DO
C	write(lout,*) 'EMISSION: ',EM_FACT(0)*C_COEFF
C
	DO NY=1,NQ(KR)
	   P=PHI(KR,NY)*C_COEFF 
	   EMISS(KR,0,NY)=EM_FACT(0)*P
	   EMISS(KR,1,NY)=EM_FACT(1)*P
	   EMISS(KR,2,NY)=EM_FACT(2)*P
	   P=PHIP(KR,NY)*C_COEFF 
	   EMISS(KR,3,NY)=-EM_FACT(3)*P
	   EMISS(KR,4,NY)=-EM_FACT(4)*P
	END DO
C  	Diagnostic output
C	WRITE(*,*) EMISS(1,0,1),EMISS(1,1,1),EMISS(1,2,1),EMISS(1,3,1),EMISS(1,4,1)
C	WRITE(*,*) EM_FACT(0),EM_FACT(1),EM_FACT(2),EM_FACT(3),EM_FACT(4)
C	WRITE(*,*) D_COEFF,SIGMA,T0(0,2),T0(1,2),T0(2,2),T0(3,1),FLARMOR,GEFF,E_COEFF
C	WRITE(*,*) SIGMA
	RETURN
	END

	SUBROUTINE KCOR(KR,EI,EQ,EU)
	INCLUDE 'PREC'
	INCLUDE 'PARAM'
	INCLUDE 'CSE'
	include 'CGRID'
	INCLUDE 'CSLINE'
	INCLUDE 'CCONST'
	INCLUDE 'CATMOS'
	INCLUDE 'CATOM'
	INCLUDE 'CORON'
	INCLUDE 'SGNM'
	INCLUDE 'CINTS'
	INCLUDE 'CLU'
C
C NEW ROUTINE OCTOBER 2019 PGJ 
C
C RETURNS K-CORONAL EMISSION COEFFICIENTS 
C
	RRR=SQRT(GX*GX+GY*GY+GZ*GZ)
	ROTANG=ATAN(GZ,GY)
	SGG=1./RRR
	CGG=SQRT(1.-SGG*SGG)
	QQ=ZERO
C	write(*,*) "RRR, ROTANG, SGG, CGG, QQ"
C	write(*,*) RRR, ROTANG, SGG, CGG, QQ
C      
C   EQUATIONS (11) AND (12) OF VANDEHULST 1950
C
        TAPB = (1.-QQ)/(1.-QQ/3.) * 2.*(1.-CGG) + 
     *           QQ/(1.-QQ/3.)      * (1.-CGG*CGG/SGG * LOG((1.+SGG)/CGG) )

	TAMB = (1.-QQ)/(1.-QQ/3.) *  2./3.*(1.-CGG**3) + 
     *       QQ/(1.-QQ/3.)       *  .25* (1.+SGG*SGG - CGG**4/SGG * 
     *                                 LOG((1.+SGG)/CGG) )

	AAA=(TAPB+TAMB)/2.
	BBB=(TAPB-TAMB)/2.
C   
C EQUATION (13) VDH 1950
C
	SIGMAT=6.6524587E-25   
	FF=CC*1.E8/ALAMB(KR)
	TR=5900.
	BETA=HH*FF/BK/TR
	B0= 2.*HH*FF * (FF/CC)**2. / (EXP(BETA) -1.)
C
C       WE NEED INTENSITY PER ANGSTROM FOR OUTPUT

C
	WWW=ALAMB(KR)
	B0 = B0 * CC*1.E8/WWW/WWW
C      
	ST=GY/SQRT(GY*GY+GX*GX)   
	CT=SQRT(1.-ST*ST)
C
C EQUATIONS (14) AND (15) VDH WHICH GIVE (16) AND (17) INTEGRANDS
C
	TTK=3./8. * B0 * SIGMAT * AAA * ED
	RRK=3./8. * B0 * SIGMAT * (AAA * CT*CT + BBB *ST*ST) *ED
	EI=TTK+RRK
	EQQ=TTK-RRK
	EUU=0.
C	WRITE(*,*)"B0,TTK,RRK, ED", B0,TTK,RRK, ED

C
C POS ROTATION  
C
	ROTANG=-2*ROTANG
	EQ=EQQ*COS(ROTANG) - EUU*SIN(ROTANG)
	EU=EQQ*SIN(ROTANG) + EUU*COS(ROTANG)

	RETURN
	END