MODULE CRYSTAL_DATA

	INTEGER, PARAMETER :: Ag3AsS3_POS     = 0
	INTEGER, PARAMETER :: AgGaS2_POS      = 1
	INTEGER, PARAMETER :: AgGaSe2_1_POS	  = 2
	INTEGER, PARAMETER :: AgGaSe2_2_POS   = 3
	INTEGER, PARAMETER :: BBO_1_POS       = 4
	INTEGER, PARAMETER :: BBO_2_POS       = 5
	INTEGER, PARAMETER :: BIBO_POS        = 6
	INTEGER, PARAMETER :: GaSe_POS        = 7
	INTEGER, PARAMETER :: K3Li2Nb5O15_POS = 8
	INTEGER, PARAMETER :: KDStarP_1_POS   = 9
	INTEGER, PARAMETER :: KDStarP_2_POS   = 10
	INTEGER, PARAMETER :: KDP_1_POS       = 11
	INTEGER, PARAMETER :: KDP_2_POS       = 12
	INTEGER, PARAMETER :: KNbO3_POS       = 13
	INTEGER, PARAMETER :: KTA_POS         = 14
	INTEGER, PARAMETER :: KTP_POS         = 15
	INTEGER, PARAMETER :: LBO_POS      = 16
	INTEGER, PARAMETER :: LiIO3_7_POS     = 17
	INTEGER, PARAMETER :: LiNbO3_POS      = 18
	INTEGER, PARAMETER :: MagnesiumBariumFlouride_POS = 19
	INTEGER, PARAMETER :: Tl3AsSe3_POS    = 20
	INTEGER, PARAMETER :: ZnGeP2_POS      = 21
	

	INTEGER, PARAMETER :: NumberOfCrystals = 22
	INTEGER CurrentCrystal

	CHARACTER*(50) CrystalName

	INTEGER CrystalType
	INTEGER, PARAMETER :: PositiveUniaxial = 1
	INTEGER, PARAMETER :: NegativeUniaxial = 2
	INTEGER, PARAMETER :: Biaxial = 3

	INTEGER CrystalClass
	INTEGER, PARAMETER :: Unknown    = 0
	INTEGER, PARAMETER :: Class_1    = 1
	INTEGER, PARAMETER :: Class_2    = 2
	INTEGER, PARAMETER :: Class_m    = 3
	INTEGER, PARAMETER :: Class_222  = 4
	INTEGER, PARAMETER :: Class_mm2  = 5
	INTEGER, PARAMETER :: Class_4    = 6
	INTEGER, PARAMETER :: Class_4_   = 7
	INTEGER, PARAMETER :: Class_4mm  = 8
	INTEGER, PARAMETER :: Class_4_2m = 9
	INTEGER, PARAMETER :: Class_3    = 10
	INTEGER, PARAMETER :: Class_32   = 11
	INTEGER, PARAMETER :: Class_3m   = 12
	INTEGER, PARAMETER :: Class_6    = 13
	INTEGER, PARAMETER :: Class_6_   = 14
	INTEGER, PARAMETER :: Class_6mm  = 15
	INTEGER, PARAMETER :: Class_6_2m = 16
	INTEGER, PARAMETER :: Class_23   = 17
	INTEGER, PARAMETER :: Class_4_3m = 18

	INTEGER XType, YType, ZType
	INTEGER, PARAMETER :: Sellmeir    = 1
	INTEGER, PARAMETER :: PowerSeries = 2
	
	REAL*8 Ax(13), Ay(13), Az(13)
	REAL*8 Abs_Min, Abs_Max
	
	REAL*8 D(3,6)

END MODULE CRYSTAL_DATA

!***************************************************************************
!
!	   FUNCTION WAVELENGTHS_VALID(CRYSTAL, Lambda_Pump, Lambda_Sgnl)
!
!	   Tests to see if the index can be calculated for all the indicated wavelengths
!
!***************************************************************************

FUNCTION WAVELENGTHS_VALID()
	USE PM_DATA; USE CRYSTAL_DATA
	IMPLICIT NONE

	LOGICAL WAVELENGTHS_VALID

	Lambda_Idlr = 1/(1/Lambda_Pump - 1/Lambda_Sgnl)
	
	IF (Lambda_Pump .GE. ABS_MIN .And. &
		Lambda_Pump .LE. ABS_MAX .And. &
		Lambda_Sgnl .GE. ABS_MIN .And. &
		Lambda_Sgnl .LE. ABS_MAX .And. &		
		Lambda_Idlr .GE. ABS_MIN .And. &
		Lambda_Idlr .LE. ABS_MAX ) THEN       
		
		WAVELENGTHS_VALID = .True.
	ELSE
		WAVELENGTHS_VALID = .False.
	ENDIF

END FUNCTION WAVELENGTHS_VALID


!************************************************************************************
!
!       SUBROUTINE SET_CRYSTAL(Crystal)
!
!		Sets D_Eff, the index coefficients, the wavelengths where they are valid,
!		and general info about the crystal
!
!************************************************************************************


SUBROUTINE SET_CRYSTAL(Crystal)
	USE CRYSTAL_DATA

	IMPLICIT NONE

!	Input

	INTEGER, INTENT(IN) :: Crystal

!	Local Variables

	INTEGER I, J

!	All the D coefficients used in this part are in 10^-12 V^-1.m 
	DO I = 1,3
		DO J = 1,6
			D(I,J) = 0.0d0 
		ENDDO
	ENDDO

	SELECT CASE (Crystal)
	CASE (Ag3AsS3_POS)
		CrystalName = 'Ag3AsS3 (Proustite)'
		CrystalType = NegativeUniaxial
		CrystalClass = Class_3m
		Abs_Min = 0.5d0
		Abs_Max = 10.0d0
!		D_Reference = 'K. F. Hulme, P.Jones, P. H. Davies, M. V. Hobden: Appl. Phys. Lett. 10, 133 (1967); 
!						D. S. Chemla, P.I. Kupecek, C. A. Schwartz: Opt. Commun. 7,225 (1973)'
		D(3,1) = 11.3d0 
		D(2,2) = 18.0d0
	CASE (AgGaS2_POS)
		CrystalName = 'AgGaS2'
		CrystalType = NegativeUniaxial
		CrystalClass = Class_4_2m
		Abs_Min = 0.5d0
		Abs_Max = 13.0d0
!		D_Reference = 'B. F. Levine, C. G. Bethea: Appl. Phys. Lett. 20, 272 (1972)'
!		D value measured at 10.6µm
		D(3,6) = 13.4d0
	CASE (AgGaSe2_1_POS)
		CrystalName = 'AgGaSe2 Ref 1'
		CrystalType = NegativeUniaxial
		CrystalClass = Class_4_2m
		Abs_Min = 0.8d0
		Abs_Max = 17.0d0
!		D_Reference = 'B. F. Levine, C. G. Bethea: Appl. Phys. Lett. 20, 272 (1972); 
!					L. D. Boyd, H. M. Kasper, J. H. McFee, F. G. Stortz: IEEE J. QE-8, 900 (1972)'
!		D value @ 10.6µm 
		D(3,6) = 33d0
	CASE (AgGaSe2_2_POS)
		CrystalName = 'AgGaSe2 Ref 2'
		CrystalType = NegativeUniaxial
		CrystalClass = Class_4_2m
		Abs_Min = 0.8d0
		Abs_Max = 17.0d0
!		D_Reference = 'B. F. Levine, C. G. Bethea: Appl. Phys. Lett. 20, 272 (1972); 
!		L. D. Boyd, H. M. Kasper, J. H. McFee, F. G. Stortz: IEEE J. QE-8, 900 (1972)'
!		D value @ 10.6µm 
		D(3,6) = 33d0
	CASE (BBO_1_POS)
		CrystalName = 'BBO Ref 1'
		CrystalType = NegativeUniaxial
		CrystalClass = Class_3m
		Abs_Min = 0.2d0
		Abs_Max = 2.1d0
		!D_Reference = '. Chen, B. Wu, A. Jiang, G. You: Sci. Sin. B 28,235 (1985)'
		D(3,1) = 0.12d0
		D(2,2) = 1.78d0
	CASE (BBO_2_POS)
		CrystalName = 'BBO Ref 2'
		CrystalType = NegativeUniaxial
		CrystalClass = Class_3m
		Abs_Min = 0.3d0
		Abs_Max = 5.0d0
		!D_Reference = '. Chen, B. Wu, A. Jiang, G. You: Sci. Sin. B 28,235 (1985)'
		D(3,1) = 0.12d0
		D(2,2) = 1.78d0
	CASE (BIBO_POS)
		CrystalName = 'BIBO'
		CrystalType = Biaxial
		CrystalClass = Unknown;
		Abs_Min = 0.286d0
		Abs_Max = 2.5d0
		D(1,2) = 2.3d0
		D(1,4) = 2.3d0
		D(2,5) = 2.4d0
		D(1,1) = 2.53d0
		D(1,3) = 1.3d0
		D(2,6) = 2.8d0
	CASE (GaSe_POS)
		CrystalName = 'GaSe'
		CrystalType = NegativeUniaxial
		CrystalClass = Class_6_2m
		Abs_Min = 0.65d0
		Abs_Max = 18.0d0
!		D_Reference = 'B. F. Levine, C. G. Bethea: Appl. Phys. Lett. 20, 272 (1972); 
!			L. D. Boyd, E. Buechler, F. G. Stortz: Appl. Phys. Lett. 18, 301 (1971)'
		D(2,2) = 54.4d0
	CASE (K3Li2Nb5O15_POS)
		CrystalName = 'K3Li2Nb5O15'
		CrystalType = NegativeUniaxial
		CrystalClass = Class_4mm
		Abs_Min = 0.4d0
		Abs_Max = 5.0d0
!		D_Reference = 'S. Singh: "Nonlinear Optical Materials", in Handbook of Lasers with Selected Data on
!			Optical Technology, ed. by R. J. Pressley (Chemical Rubber Co., Cleveland, OH 1971)'
		D(1,5) = 6.2d0
		D(3,1) = 7.0d0
		D(3,3) = 12.7d0
!		(d33 never appear in Deff)
!		For this Crystal, we do not have the Kleinman symmetry. 
!		We have calculate the Deff by using an e-ray for the pump
	CASE (KDStarP_1_POS)
		CrystalName = 'KD*P Ref 1'
		CrystalType = NegativeUniaxial
		CrystalClass = Class_4_2m
		Abs_Min = 0.2d0
		Abs_Max = 1.8d0
!		D_Reference = 'B. F. Levine, C. G. Bethea: Appl. Phys. Lett. 20, 272 (1972); 
!			S. Kielich: Opto-electronics 2, 125 (1970)'
!		D value measured at 1.06µm
		D(3,6) = 0.402d0
	CASE (KDStarP_2_POS)
		CrystalName = 'KD*P Ref 2'
		CrystalType = NegativeUniaxial
		CrystalClass = Class_4_2m
		Abs_Min = 0.2d0
		Abs_Max = 1.8d0
!		D_Reference = 'B. F. Levine, C. G. Bethea: Appl. Phys. Lett. 20, 272 (1972); 
!			S. Kielich: Opto-electronics 2, 125 (1970)'
!		D value measured @ 1.06µm
		D(3,6) = 0.402d0
	CASE (KDP_1_POS)
		CrystalName = 'KDP Ref 1'
		CrystalType = NegativeUniaxial
		CrystalClass = Class_4_2m
		Abs_Min = 0.2d0
		Abs_Max = 1.5d0
!		D_Reference = 'B. F. Levine, C. G. Bethea: Appl. Phys. Lett. 20, 272 (1972)'
!		D value @ 1.06µm
		D(1,4) = 0.435d0
		D(2,5) = 0.435d0
		D(3,6) = 0.435d0
	CASE (KDP_2_POS)
		CrystalName = 'KDP Ref 2'
		CrystalType = NegativeUniaxial
		CrystalClass = Class_4_2m
		Abs_Min = 0.2d0
		Abs_Max = 2.5d0
!		D_Reference = 'B. F. Levine, C. G. Bethea: Appl. Phys. Lett. 20, 272 (1972)'
!		D value @ 1.06µm
		D(1,4) = 0.435d0
		D(2,5) = 0.435d0
		D(3,6) = 0.435d0
	CASE (KNbO3_POS)
		CrystalName = 'KNbO3'
		CrystalType = Biaxial
		CrystalClass = Unknown
		Abs_Min = 0.4d0
		Abs_Max = 4.5d0
	CASE (KTA_POS)
		CrystalName = 'KTA'
		CrystalType = Biaxial
		CrystalClass = Unknown
		Abs_Min = 0.35d0
		Abs_Max = 4.0d0
	CASE (KTP_POS)
		CrystalName = 'KTP'
		CrystalType = Biaxial
		CrystalClass = Unknown
		Abs_Min = 0.35d0
		Abs_Max = 4.50d0
	
	CASE (LBO_POS)
		CrystalName = 'LBO'
		CrystalType = Biaxial
		CrystalClass = Unknown
		Abs_Min = 0.365d0
		Abs_Max = 1.1d0
	
	CASE (LiIO3_7_POS)
		CrystalName = 'LiIO3 Ref 7'
		CrystalType = NegativeUniaxial
		CrystalClass = Class_6
		Abs_Min = 0.3d0
		Abs_Max = 5.0d0
!		D_Reference = 'B. F. Levine, C. G. Bethea: Appl. Phys. Lett. 20, 272 (1972)'
!		D values measured at 1.06µm
		D(1,5) = -5.53d0
		D(2,4) = -5.53d0
		D(3,1) = -5.53d0
		D(3,2) = -5.53d0
		D(3,3) = -5.53d0
	CASE (LiNbO3_POS)
		CrystalName = 'LiNbO3'
		CrystalType = NegativeUniaxial
		CrystalClass = Class_3m
		Abs_Min = 0.4d0
		Abs_Max = 3.4d0
!		D_Reference = 'B. F. Levine, C. G. Bethea: Appl. Phys. Lett. 20, 272 (1972)'
!		D values measured at 1.06µm
		D(1,5) = -5.44d0
		D(2,2) = 2.76d0
		D(3,1) = -5.44d0
	CASE (MagnesiumBariumFlouride_POS)
		CrystalName = 'Magnesium-Barium Flouride'
		CrystalType = Biaxial
		CrystalClass = Unknown
		Abs_Min = 0.3d0
		Abs_Max = 5.0d0
	CASE (Tl3AsSe3_POS)
		CrystalName = 'Tl3AsSe3'
		CrystalType = NegativeUniaxial
		CrystalClass = Class_3m
		Abs_Min = 1.26d0
		Abs_Max = 16.0d0
!		D_Reference = 'Feichtner, J.D., G.W. Roland, Appl.Opt., 11, 993 (1972)'
!		Conversion of the units of nonlinear coefficients from the SI system to CGS system can be 
!		performed if we bear in mind that 1 m/V=2.39*10^3 esu
		D(3,1) = 12.5d0
		D(2,2) = 13.4d0
	CASE (ZnGeP2_POS)
		CrystalName = 'ZnGeP2'
		CrystalType = PositiveUniaxial
		CrystalClass = Class_4_2m
		Abs_Min = 0.74d0
		Abs_Max = 12.0d0
!		D_Reference = 'B. F. Levine, C. G. Bethea: Appl. Phys. Lett. 20, 272 (1972)'
!		D value measure at 10.6µm
		D(3,6) = 75.4d0
	CASE DEFAULT
		CrystalName = 'No Crystal Selected'
		CrystalType = Unknown
		CrystalClass = Unknown
		Abs_Min = 0.1
		Abs_Max = 0.2
	END SELECT
END SUBROUTINE SET_CRYSTAL

SUBROUTINE INDEX(Crystal, Lambda, T, Nx, Ny, Nz)
	USE CRYSTAL_DATA
	IMPLICIT NONE

!	Arguments

	INTEGER, INTENT(IN) :: Crystal
	REAL*8, INTENT(IN) :: Lambda, T
	REAL*8, INTENT(OUT) :: Nx, Ny, Nz

!	Local Variables
	REAL*8 T1, T2, T3
	REAL*8 X1, X2, X3, Y1, Y2, Y3, Z1, Z2, Z3

!	Begin Subroutine

	SELECT CASE (Crystal)
	CASE (Ag3AsS3_POS)
!		Ag3AsS3 ( proustite )
!		G. AIROLDI, "Nonlinear Materials for Tunable Coherent Light 
!		Sources in the Middle Infra-Red," Rivista Del Nuovo Cimento,
!		Vol 6, N 3, p 295-319 (1976)
		
		Nx = DSQRT( 9.220 + .4454 / ( Lambda**2 - 0.1264 )	&
					+ 1733. / (Lambda**2 - 1000.) )
		Ny = Nx
		Nz = DSQRT( 7.007 + .3230 / ( Lambda**2 - 0.1192 )	&
					+ 660. / ( Lambda**2 - 1000.) )

	CASE (AgGaS2_POS)
!		AgGaS2
!		G. C. Bahr, Appl. Opt. 15, 305 (1976)

		Nx = DSQRT( 3.6280 + 2.1686*Lambda**2/(Lambda**2 - 0.1003)	&
					+ 2.1753*Lambda**2/(Lambda**2 - 950 ) )
		Ny = Nx
		Nz = DSQRT( 4.0172+ 1.5274*Lambda**2/(Lambda**2 - 0.1310)	&
					+ 2.1699*Lambda**2/(Lambda**2 - 950) )

	CASE (AgGaSe2_1_POS)
!		AgGaSe2 Ref 1
!		Kildal, H., J. Mikkelsen, Opt. Commun., 9, 315, (1973)
!		~0.8 - 17 (35% Refl. loss)
!		1 - 13.5

		Nx = DSQRT( 3.9362 + 2.9113/(1 - (0.38821/Lambda)**2) +	 &
             1.7954/(1 - (40/Lambda)**2) )
		Ny = Nx
		Nz = DSQRT( 3.3132 + 3.3616/(1 - (0.38201/Lambda)**2) +	 &
             1.7677/(1 - (40/Lambda)**2) )

	CASE (AgGaSe2_2_POS)
!		AgGaSe2 Ref 2
!		Bhar, G. C., Appl. Opt., 15, 305 (1976)
!		~0.8 - 17 (35% Refl. loss)
!		0.7 - 13.5

		Nx = DSQRT( 4.6453 + 2.2057/(1 - (0.43347/Lambda)**2) +	  &
			1.8377/(1-(40/Lambda)**2) )
		Ny = Nx
		Nz = DSQRT( 5.2912 + 1.3970/(1 - (0.53339/Lambda)**2) +	  &
			1.9282/(1 - (40/Lambda)**2) )						  

	CASE (BBO_1_POS)
!		BBO Ref 1
!		Katz, K., IEEE J. Quant. Electron., 22, 1013-4 (1986)
!		~0.2 - 2.1
!		0.3 - 1.5 ?

		Nx = DSQRT( 2.7359 + 0.01878/(Lambda**2 - 0.01822) -	   &
             0.01354*Lambda**2 )
		Ny = Nx
		Nz = DSQRT( 2.3753 + 0.01224/(Lambda**2 - 0.01667) -	   &
             0.01516*Lambda**2 )

	CASE (BBO_2_POS)
!		BBO Ref 2
!		Graham, A. Zalkin, J. Appl. Phys., 62 (1987)
!		0.2-2.8
!		0.4 - 1.0

		IF ( Lambda .LE. 12.0 ) THEN
			Nx = DSQRT( 2.7405 + 0.0184/(Lambda**2 - 0.0179) -	 &
				0.0155*Lambda**2 )
			Ny = Nx
			Nz = DSQRT( 2.3730 + 0.0128/(Lambda**2 - 0.0156) -	 &
				0.0044*Lambda**2 )
		ELSE
			Nx = 0.2
			Ny = Nx
			Nz = 2.8
		ENDIF

	CASE (BIBO_POS)

		Nx = DSQRT(3.0740d0 + 0.0323d0/(Lambda**2-0.0316d0) - 0.01337d0*Lambda**2 )
		Ny = DSQRT(3.1685d0 + 0.0373d0/(Lambda**2-0.0346d0) - 0.01750d0*Lambda**2 )
		Nz = DSQRT(3.6545d0 + 0.0511d0/(Lambda**2-0.0371d0) - 0.0226d0*Lambda**2  )

	CASE (GaSe_POS)
!		GaSe
!		Abdullaev, G.B., L.A. Kulevskii, A.M. Prokharov, A.D. Savel'ev,
!		E. Yu. Salaev, V.V. Smirnov, JETP Lett., 16, 90 (1972)
!		~0.65 - 18 (40% Relf. loss)
!		0.63 - 10.6

		Nx = DSQRT( -0.05466/Lambda**4 + 0.48605/Lambda**2 +	 &
             7.8902 - (0.000824)*Lambda**2 -				 &
             (0.00000273)*Lambda**4 )
		Ny = Nx
		Nz = DSQRT( 6.0476 + 0.3423/(Lambda**2 - 0.16491) -	  &
             (0.001042)*Lambda**2 )

	CASE (K3Li2Nb5O15_POS)
!		K3Li2Nb5O15
!		Nikogosyan, Handbook of Nonlinear Optical Crystals, 1991, p. 76
!		0.45 - 1?

		Nx = DSQRT( 1.0 + 3.708*Lambda**2/(Lambda**2 - 0.04601) )
		Ny = Nx
		Nz = DSQRT( 1.0 + 3.349*Lambda**2/(Lambda**2 - 0.03564) )
	CASE (KDStarP_1_POS)
!		KD*P Ref 1
!		Barnes, N.P., D.J. Gettemy, R.S. Adhav, J. Opt. Soc. Am.,
!		72, 895 (1982)
!		~0.2 - 1.8 (98% composition)
!		0.41 - 0.65

		Nx = DSQRT( 1.0 + 1.23318*Lambda**2/(Lambda**2 - 0.00931) )
		Ny = Nx
		Nz = DSQRT( 1.0 + 1.12354*Lambda**2/(Lambda**2 - 0.00828) )

	CASE (KDStarP_2_POS)
!		KD*P Ref 2
!		Phillips, R.A., J. Opt. Soc. Am., 56, 629 (1966)
!		~0.2 - 1.8 (98% composition)

		Nx = DSQRT( 1.661824 + 0.585337/(1 - (0.016017/Lambda**2)) +	  &
			0.691221/(1 - (30.0/Lambda**2)) )
		Ny = Nx
		Nz = DSQRT( 1.687522 + 0.447488/(1 - (0.017039/Lambda**2)) +	  &
			0.596216/(1 - (30.0/Lambda**2)) )

	CASE (KDP_1_POS)
!		KDP Ref 1
!		Zernike, Frits, Jr., J. Opt. Soc. Am., 54, 1215 (1964)
!		~0.2 - 1.5
!		0.21 - 1.5

		Nx = DSQRT( 0.01011279/(Lambda**2 - 0.01294238) +		&
			0.03249268/(0.0025 - 1/Lambda**2) + 2.260476 )
		Ny = Nx
		Nz = DSQRT( 0.00865325/(Lambda**2 - 0.01229326) +		&
             0.00806984/(0.0025 - 1/Lambda**2) + 2.133831 )

	CASE (KDP_2_POS)
!		KDP Ref 2
!		Barnes, N.P., D.J. Gettemy, R.S. Adhav, J. Opt. Soc. Am.,
!		72, 895 (1982)
!		0.41 - 0.65

		Nx = DSQRT( 1.0 + 1.24361*Lambda**2/(Lambda**2 - 0.00959) )
		Ny = Nx
		Nz = DSQRT( 1.0 + 1.12854*Lambda**2/(Lambda**2 - 0.00841) )

	CASE (KNbO3_POS)
!		KNbO3
!		J. Opt. Soc. Am. B/Vol.9,No.3/March 1992
!		Refractive indices of orthorhombic KNbO3
!		B. Zysset
		IF (T .LT. 75d0) THEN
			T1 = 22
			T2 = 50
			T3 = 75
		ELSEIF (T .LT. 100d0) THEN
			T1 = 50
			T2 = 75
			T3 = 100
		ELSEIF (T .LT. 140d0) THEN
			T1 = 75
			T2 = 100
			T3 = 140
		ELSE
			T1 = 100
			T2 = 140
			T3 = 180
		ENDIF
		CALL KNBO3_INDEX(Lambda, T1, X1, Y1, Z1)
		CALL KNBO3_INDEX(Lambda, T2, X2, Y2, Z2)
		CALL KNBO3_INDEX(Lambda, T3, X3, Y3, Z3)

!		Quadratic Interpolation / Extrapolation
		Nx = X1*(T-T2)*(T-T3)/(T1-T2)/(T1-T3) + &
			 X2*(T-T1)*(T-T3)/(T2-T1)/(T2-T3) + &
			 X3*(T-T1)*(T-T2)/(T3-T1)/(T3-T2) 
		Ny = Y1*(T-T2)*(T-T3)/(T1-T2)/(T1-T3) + &
			 Y2*(T-T1)*(T-T3)/(T2-T1)/(T2-T3) + &
			 Y3*(T-T1)*(T-T2)/(T3-T1)/(T3-T2) 
		Nz = Z1*(T-T2)*(T-T3)/(T1-T2)/(T1-T3) + &
			 Z2*(T-T1)*(T-T3)/(T2-T1)/(T2-T3) + &
			 Z3*(T-T1)*(T-T2)/(T3-T1)/(T3-T2) 

	CASE (KTA_POS)
!		KTA
!		D. L. Fenimore, K. L. Schepler, U. B. Ramabadran,
!		S. R. McPherson, "Infrared corrected Sellmeier coefficients
!		for potassium titaNyl arsenate," J. Opt. Soc. B, 12 (5),
!		794-796 (May 1995)

		Nx = DSQRT( 1.90713 + 1.23522/(1 - (0.19692/Lambda)**2) - (0.01025*Lambda**2))
		Ny = DSQRT( 2.15912 + 1.00099/(1 - (0.21844/Lambda)**2) - (0.01096*Lambda**2))
		Nz = DSQRT( 2.14786 + 1.29559/(1 - (0.22719/Lambda)**2) - (0.01436*Lambda**2))

	CASE (KTP_POS)
!		KTP
!		H. Vanherzeele, J. D. Bierlein, F. C. Zumsteg, Appl. Opt., 27,
!		3314 (1988)

		Nx = DSQRT( 2.1146 + 0.89188/(1 - (0.20861/Lambda)**2) - (0.01320*Lambda**2))
		Ny = DSQRT( 2.1518 + 0.87862/(1 - (0.21801/Lambda)**2) - (0.01327*Lambda**2))
		Nz = DSQRT( 2.3136 + 1.00012/(1 - (0.23831/Lambda)**2) - (0.01679*Lambda**2))

	CASE (LBO_POS)
!		LBO
!		Velsko et al, "Phase-matched harmonic generation in Lithum Triborate (LBO)"
!		IEEE J. Quant. Eletron. 27, 2182-2192(1991)
!		0.365 - > 1.1 microns


		Nx = DSQRT( 2.4543d0 + .011413d0/(Lambda**2-0.0094981) - .0138900d0*Lambda**2 )
		Ny = DSQRT( 2.5382d0 + .012830d0/(Lambda**2-0.01387) - .017034d0*Lambda**2 )
		Nz = DSQRT( 2.5854d0 + .013065d0/(Lambda**2-0.011617) - .018146d0*Lambda**2 )



	CASE (LiIO3_7_POS)
!		LiIO3
!		Phillips, R.A., J. Opt. Soc. Am., 56, 629 (1966)
!		0.3 - > 5.0
!		0.35 - 5.0

		Nx = DSQRT( 2.03132d0 + 1.37623d0/(1.0d0 - (0.0350832d0/ Lambda**2)) + 1.06745d0/(1.0d0 - (169.0d0/Lambda**2)) )
		Ny = Nx
		Nz = DSQRT( 1.83086d0 + 1.08807d0/(1.0d0 - (0.031381d0 / Lambda**2)) + 0.554582d0/(1.0d0 - (158.76d0/Lambda**2)) )

	CASE (LiNbO3_POS)
!		LiNbO3
!		Smith, D.S., H.D. Riccius, R.P. Edwin, Opt. Commun., 17,
!		332, (1976)
!		Probably to 3.4
!		0.4 - 3.39

		Nx = DSQRT( 4.9048 - 0.11768/(0.04750 - Lambda**2) -	  &
			0.027169*Lambda**2 )
		Ny = Nx
		Nz = DSQRT( 4.5820 - 0.099169/(0.044432 - Lambda**2) -  &
			0.021950*Lambda**2 )

	CASE (MagnesiumBariumFlouride_POS)
!		Magnesium-Barium Flouride
!		V. G. Dmitriev, G. G. Gurzadyan, D. N. Nikogosyan,
!		"Handbook of Nonlinear Crystals," (Springer-Verlag,
!		Berlin) 1991, pp.112-113
		Nx = DSQRT( 2.077 + 0.0076/(Lambda**2 - 0.0079) )
		Ny = DSQRT( 2.1238 + 0.0086/(Lambda**2) )
		Nz = DSQRT( 2.1462 + 0.00736/(Lambda**2 - 0.009) )

	CASE (Tl3AsSe3_POS)
!		Tl3AsSe3
!		Feichtner, J.D., G.W. Roland, Appl.Opt., 11, 993 (1972)
!		1.26 - 16 (45% not corrected for refl.)
!		1.5 - 5.3

		Nx = DSQRT( 1.0 + 9.977*Lambda**2/(Lambda**2 - (0.435)**2) +	&
			0.067*Lambda**2/(Lambda**2 - (20)**2) )
		Ny = Nx
		Nz = DSQRT( 1.0 + 8.783*Lambda**2/(Lambda**2 - (0.435)**2) +  &
			0.051*Lambda**2/(Lambda**2 - (20)**2) )

	CASE (ZnGeP2_POS)
!		ZnGeP2
!		Dmitriev p. 85
!		0.74 - 12

		Nx = DSQRT( 4.4733 + 5.26575*Lambda**2/(Lambda**2 - 0.13381) +   &
			1.49085*Lambda**2/(Lambda**2 - 662.55 ) )
		Ny = Nx
		Nz = DSQRT( 4.63318 + 5.34215*Lambda**2/(Lambda**2 - 0.142555) +	&
			1.45785*Lambda**2/(Lambda**2 - 662.55) )

	CASE DEFAULT
		Nx = 1d0
		Ny = 1d0
		Nz = 1d0
	END SELECT

END SUBROUTINE INDEX


REAL*8 FUNCTION SELLMEIR1(L, L1, L2, S1, S2, D)

	REAL*8, INTENT(IN) :: L, L1, L2, S1, S2, D

	REAL*8 Arg

	Arg = 1d0 + S1 * L**2 * L1**2 / (L**2-L1**2) &
						  + S2 * L**2 * L2**2 / (L**2-L2**2) &
						  - D * L**2
	SELLMEIR1 = DSQRT(Arg)
END FUNCTION SELLMEIR1

SUBROUTINE KNBO3_INDEX(Lambda, T, Nx, Ny, Nz)
	IMPLICIT NONE
	
!	Arguments
	
	REAL*8, INTENT(IN) :: Lambda, T
	REAL*8, INTENT(OUT) :: Nx, Ny, Nz

!	Local Variables

	REAL*8 D1, L1, L2, S1, S2, SELLMEIR1

!	Begin Subroutine

	IF (T .EQ. 22) THEN
		D1 = 1.943289D-2
		L1 = 2.55227711D-01
		L2 = 1.19714173D-01
		S1 = 1.609170D+01
		S2 = 1.654431D+02
		Nx = SELLMEIR1(Lambda,L1,L2,S1,S2,D1)

		D1 = 2.517432D-2
		L1 = 2.58157325D-01
		L2 = 1.29092115D-01
		S1 = 2.005519D+01
		S2 = 1.498408D+02
		Ny = SELLMEIR1(Lambda,L1,L2,S1,S2,D1)
							
		D1 = 2.845018D-2
		L1 = 2.72745560D-01
		L2 = 1.37003854D-01
		S1 = 1.937347D+01
		S2 = 1.354992D+02
		Nz = SELLMEIR1(Lambda,L1,L2,S1,S2,D1)
	ELSEIF (T .EQ. 50) THEN
		D1 = 1.945428D-2
		L1 = 2.52752522D-01
		L2 = 1.11742242D-01
		S1 = 1.824447D+01
		S2 = 1.811303D+02
		Nx = SELLMEIR1(Lambda,L1,L2,S1,S2,D1)

		D1 = 2.514857D-2
		L1 = 2.58664133D-01
		L2 = 1.27592173D-01
		S1 = 2.027475D+01
		S2 = 1.523380D+02
		Ny = SELLMEIR1(Lambda,L1,L2,S1,S2,D1)

		D1 = 2.827127D-2
		L1 = 2.71398770D-01
		L2 = 1.30593386D-01
		S1 = 2.066250D+01
		S2 = 1.441607D+02
		Nz = SELLMEIR1(Lambda,L1,L2,S1,S2,D1)
	ELSEIF (T .EQ. 75) THEN
		D1 = 1.952880D-2
		L1 = 2.54675819D-01
		L2 = 1.13972981D-01
		S1 = 1.756954D+01
		S2 = 1.766712D+02
		Nx = SELLMEIR1(Lambda,L1,L2,S1,S2,D1)

		D1 = 2.514893D-2
		L1 = 2.60772236D-01
		L2 = 1.29504861D-01
		S1 = 1.933691D+01
		S2 = 1.504971D+02
		Ny = SELLMEIR1(Lambda,L1,L2,S1,S2,D1)

		D1 = 2.807510D-2
		L1 = 2.72530570D-01
		L2 = 1.31949811D-01
		S1 = 2.006891D+01
		S2 = 1.427868D+02
		Nz = SELLMEIR1(Lambda,L1,L2,S1,S2,D1)
	ELSEIF (T .EQ. 100) THEN
		D1 = 1.968455D-2
		L1 = 2.62219697D-01
		L2 = 1.25100015D-01
		S1 = 1.401940D+01
		S2 = 1.584156D+02
		Nx = SELLMEIR1(Lambda,L1,L2,S1,S2,D1)

		D1 = 2.515860D-2
		L1 = 2.65628695D-01
		L2 = 1.37103071D-01
		S1 = 1.660832D+01
		S2 = 1.420508D+02
		Ny = SELLMEIR1(Lambda,L1,L2,S1,S2,D1)

		D1 = 2.796330D-2
		L1 = 2.81955899D-01
		L2 = 1.51431159D-01
		S1 = 1.437713D+01
		S2 = 1.233481D+02
		Nz = SELLMEIR1(Lambda,L1,L2,S1,S2,D1)
	ELSEIF (T .EQ. 140) THEN
		D1 = 1.982648D-2
		L1 = 2.58712656D-01
		L2 = 1.18491161D-01
		S1 = 1.648284D+01
		S2 = 1.677579D+02
		Nx = SELLMEIR1(Lambda,L1,L2,S1,S2,D1)

		D1 = 2.508202D-2
		L1 = 2.67779444D-01
		L2 = 1.38395441D-01
		S1 = 1.604692D+01
		S2 = 1.407942D+02
		Ny = SELLMEIR1(Lambda,L1,L2,S1,S2,D1)

		D1 = 2.771253D-2
		L1 = 2.74864037D-01
		L2 = 1.35282050D-01
		S1 = 1.887445D+01
		S2 = 1.388675D+02
		Nz = SELLMEIR1(Lambda,L1,L2,S1,S2,D1)
	ELSEIF (T .EQ. 180) THEN
		D1 = 2.004523D-2
		L1 = 2.52773288D-01
		L2 = 9.88011467D-02
		S1 = 2.147866D+01
		S2 = 2.156493D+02
		Nx = SELLMEIR1(Lambda,L1,L2,S1,S2,D1)

		D1 = 2.511282D-2
		L1 = 2.73698130D-01
		L2 = 1.43185621D-01
		S1 = 1.389662D+01
		S2 = 1.372157D+02
		Ny = SELLMEIR1(Lambda,L1,L2,S1,S2,D1)

		D1 = 2.762299D-2
		L1 = 2.71016125D-01
		L2 = 1.19793058D-01
		S1 = 2.212020D+01
		S2 = 1.628852D+02
		Nz = SELLMEIR1(Lambda,L1,L2,S1,S2,D1)
	ELSE
		Nx = 0d0
		Ny = 0d0
		Nz = 0d0
	ENDIF

END SUBROUTINE KNBO3_INDEX


!****************************************************************************************************
!
!		SUBROUTINE EFFICIENAy(Crystal,TypePM,Phi_Pump,Theta_Pump,Theta_Sgnl,Theta_Idlr,Phi_Sgnl,Phi_Idlr,Deff)
!
!		Calculates D_eff for given parameters
!
!****************************************************************************************************

REAL*8 FUNCTION GET_D_EFF()
	USE CRYSTAL_DATA; USE PM_DATA
	IMPLICIT NONE

!	Local Variables

	REAL*8 v,w
	REAL*8 a,b,c
	REAL*8 Deff

	IF ((CrystalType .EQ. PositiveUniaxial) .OR. (CrystalType .EQ. NegativeUniaxial)) THEN
!			For this calculation of the Deff we assume that Kleinman's symmetry condition is valid. 
!			So, it makes no difference which of the three waves is the output.
!			That's true when all frequencies lie in a region without absorption.
!			That is, if Kleinman's condition is taken into account, we find one expression for an expression
!			for an interaction between two o-rays and an e-ray and one for an interaction between two o-rays 
!			and an o-ray.

!			For uniaxial Crystals
!			The "fast" index, is the smaller of the two indices. So, we have to look if it's positive or a 
!			negative  unixial Crystal.

		v = DATAN(  (  DTAN(Phi_Pump)*DCOS(Phi_Sgnl)*DTAN(Theta_Sgnl)*DCOS(Theta_Pump) 			&
				  +DSIN(Phi_Sgnl)*DTAN(Theta_Sgnl) + DSIN(Theta_Pump)*DTAN(Phi_Pump) )		&
			  / ( -DCOS(Phi_Sgnl)*DTAN(Theta_Sgnl)*DCOS(Theta_Pump) 						&
			      +DSIN(Phi_Sgnl)*DTAN(Theta_Sgnl)*DTAN(Phi_Pump)-DSIN(Theta_Pump) ) )

		w = DATAN(  ( DTAN(Phi_Pump)*DCOS(Phi_Idlr)*DTAN(Theta_Idlr)*DCOS(Theta_Pump) +			&
				  DSIN(Phi_Idlr)*DTAN(Theta_Idlr)+DSIN(Theta_Pump)*DTAN(Phi_Pump) )			&
			  / (-DCOS(Phi_Idlr)*DTAN(Theta_Idlr)*DCOS(Theta_Pump) +						&
			      DSIN(Phi_Idlr)*DTAN(Theta_Idlr)*DTAN(Phi_Pump)-DSIN(Theta_Pump) ) )

		a = -1/(DSQRT((-1/DTAN(v)*DSIN(Theta_Sgnl)*DCOS(Phi_Sgnl)*DCOS(Theta_Pump)*DCOS(Phi_Pump)+ &
			1/DTAN(v)*DSIN(Theta_Sgnl)*DSIN(Phi_Sgnl)*DSIN(Phi_Pump)-1/DTAN(v)*DCOS(Theta_Sgnl)* &
			DSIN(Theta_Pump)*DCOS(Phi_Pump)+DSIN(Theta_Sgnl)*DCOS(Phi_Sgnl)*DCOS(Theta_Pump)*DSIN(Phi_Pump)+ &
			DSIN(Theta_Sgnl)*DSIN(Phi_Sgnl)*DCOS(Phi_Pump)+DCOS(Theta_Sgnl)*DSIN(Theta_Pump)* &
			DSIN(Phi_Pump))**2/(DSIN(Theta_Sgnl)*DCOS(Phi_Sgnl)*DSIN(Theta_Pump)-DCOS(Theta_Sgnl)* &
			DCOS(Theta_Pump))**2+1+1/(DTAN(v)**2))*DTAN(v))

		b = 1/(DSQRT((-1/DTAN(v)*DSIN(Theta_Sgnl)*DCOS(Phi_Sgnl)*DCOS(Theta_Pump)*DCOS(Phi_Pump)+ &
			1/DTAN(v)*DSIN(Theta_Sgnl)*DSIN(Phi_Sgnl)*DSIN(Phi_Pump)-1/DTAN(v)*DCOS(Theta_Sgnl)* &
			DSIN(Theta_Pump)*DCOS(Phi_Pump)+DSIN(Theta_Sgnl)*DCOS(Phi_Sgnl)*DCOS(Theta_Pump)*DSIN(Phi_Pump)+ &
			DSIN(Theta_Sgnl)*DSIN(Phi_Sgnl)*DCOS(Phi_Pump)+DCOS(Theta_Sgnl)*DSIN(Theta_Pump)* &
			DSIN(Phi_Pump))**2/(DSIN(Theta_Sgnl)*DCOS(Phi_Sgnl)*DSIN(Theta_Pump)-DCOS(Theta_Sgnl)* &
			DCOS(Theta_Pump))**2+1+1/(DTAN(v)**2)))

		c = DSQRT(1-1/(((-1/DTAN(v)*DSIN(Theta_Sgnl)*DCOS(Phi_Sgnl)*DCOS(Theta_Pump)*DCOS(Phi_Pump)+ &
			1/DTAN(v)*DSIN(Theta_Sgnl)*DSIN(Phi_Sgnl)*DSIN(Phi_Pump)-1/DTAN(v)*DCOS(Theta_Sgnl)* &
			DSIN(Theta_Pump)*DCOS(Phi_Pump)+DSIN(Theta_Sgnl)*DCOS(Phi_Sgnl)*DCOS(Theta_Pump)*DSIN(Phi_Pump)+ &
			DSIN(Theta_Sgnl)*DSIN(Phi_Sgnl)*DCOS(Phi_Pump)+DCOS(Theta_Sgnl)*DSIN(Theta_Pump)* &
			DSIN(Phi_Pump))**2/(DSIN(Theta_Sgnl)*DCOS(Phi_Sgnl)*DSIN(Theta_Pump)-DCOS(Theta_Sgnl)* &
			DCOS(Theta_Pump))**2+1+1/(DTAN(v)**2))*DTAN(v)**2)-1/((-1/DTAN(v)*DSIN(Theta_Sgnl)*DCOS(Phi_Sgnl)* &
			DCOS(Theta_Pump)*DCOS(Phi_Pump)+1/DTAN(v)*DSIN(Theta_Sgnl)*DSIN(Phi_Sgnl)*DSIN(Phi_Pump)- &
			1/DTAN(v)*DCOS(Theta_Sgnl)*DSIN(Theta_Pump)*DCOS(Phi_Pump)+DSIN(Theta_Sgnl)*DCOS(Phi_Sgnl)* &
			DCOS(Theta_Pump)*DSIN(Phi_Pump)+DSIN(Theta_Sgnl)*DSIN(Phi_Sgnl)*DCOS(Phi_Pump)+DCOS(Theta_Sgnl)* &
			DSIN(Theta_Pump)*DSIN(Phi_Pump))**2/(DSIN(Theta_Sgnl)*DCOS(Phi_Sgnl)*DSIN(Theta_Pump)- &
			DCOS(Theta_Sgnl)*DCOS(Theta_Pump))**2+1+1/(DTAN(v)**2)))

		IF ( (CrystalClass .EQ. Class_4) .OR. (CrystalClass .EQ. Class_6) ) THEN
		   
			IF ( TypePM .EQ. 1 ) THEN

				Deff = DABS( D(3,1) * DSIN(Theta_Pump) * DCOS(v-w) )
				  
			ELSE
			      
				Deff = DABS(- D(3,1) * (- c * DCOS(Theta_Pump) * DSIN(Phi_Pump+w) &
										+ a * DSIN(Theta_Pump) * DSIN(w) &
										+ b * DSIN(Theta_Pump) * DCOS(w) ) )
			END IF

		ELSE IF (CrystalClass .EQ. Class_32) THEN
	
			IF ( TypePM .EQ. 1 ) THEN

				Deff = DABS(-D(1,1) * DCOS(Theta_Pump) * (DCOS(Phi_Pump) * DSIN(v) * DSIN(w) &
							-DCOS(Phi_Pump) * DCOS(v) * DCOS(w) &
							-DSIN(Phi_Pump) * DSIN(v) * DCOS(w) &
							-DSIN(Phi_Pump) * DCOS(v) * DSIN(w) ) )
			ELSE

				Deff = DABS(-D(1,1) * a * DCOS(Theta_Pump) * (DCOS(Phi_Pump) * DSIN(w) &
							-b * DCOS(Phi_Pump) * DCOS(w) &
							-a * DSIN(Phi_Pump) * DCOS(w) &
							-b * DSIN(Phi_Pump) * DSIN(w) ) )
			END IF

		ELSE IF ( CrystalClass .EQ. Class_3) THEN

			IF ( TypePM .EQ. 1 ) THEN

				Deff = DABS( D(1,1) * DCOS(Theta_Pump) * DCOS(Phi_Pump-v-w) &
							-D(2,2) * DCOS(Theta_Pump) * DSIN(Phi_Pump-v-w) &
							+D(3,1) * DSIN(Theta_Pump) * DCOS(v-w) )
			ELSE

				Deff = DABS(-D(1,1) * a * DCOS(Theta_Pump) * DCOS(Phi_Pump) * DSIN(w) &
							+D(1,1) * b * DCOS(Theta_Pump) * DCOS(Phi_Pump) * DCOS(w) &
							-D(3,1) * c * DCOS(Theta_Pump) * DCOS(Phi_Pump) * DSIN(w) &
							+D(2,2) * a * DCOS(Theta_Pump) * DCOS(Phi_Pump) * DCOS(w) &
							+D(2,2) * b * DCOS(Theta_Pump) * DCOS(Phi_Pump) * DSIN(w) &
							+D(2,2) * a * DCOS(Theta_Pump) * DSIN(Phi_Pump) * DSIN(w) &
							-D(2,2) * b * DCOS(Theta_Pump) * DSIN(Phi_Pump) * DCOS(w) &
							-D(3,1) * c * DCOS(Theta_Pump) * DSIN(Phi_Pump) * DCOS(w) &
							+D(1,1) * a * DCOS(Theta_Pump) * DSIN(Phi_Pump) * DCOS(w) &
							+D(1,1) * b * DCOS(Theta_Pump) * DSIN(Phi_Pump) * DSIN(w) &
							+D(3,1) * a * DSIN(Theta_Pump) * DSIN(w) &
							+D(3,1) * b * DSIN(Theta_Pump) * DCOS(w) )
			END IF

		ELSE IF ( CrystalClass .EQ. Class_3m) THEN

			IF ( TypePM .EQ. 1 ) THEN
				Deff = DABS(  D(3,1) * DSIN(Theta_Pump) * DCOS(v-w) &
							+ D(2,2) * DCOS(Theta_Pump) * DSIN(v+w-Phi_Pump) )
			ELSE
				Deff = DABS(- D(3,1) * c * DCOS(Theta_Pump) * DCOS(Phi_Pump) * DSIN(w) &
							+ D(2,2) * a * DCOS(Theta_Pump) * DCOS(Phi_Pump) * DCOS(w) &
							+ D(2,2) * b * DCOS(Theta_Pump) * DCOS(Phi_Pump) * DSIN(w) &
							+ D(2,2) * a * DCOS(Theta_Pump) * DSIN(Phi_Pump) * DSIN(w) &
							- D(2,2) * b * DCOS(Theta_Pump) * DSIN(Phi_Pump) * DCOS(w) &
							- D(3,1) * c * DCOS(Theta_Pump) * DSIN(Phi_Pump) * DCOS(w) &
							+ D(3,1) * a * DSIN(Theta_Pump) * DSIN(w) &
							+ D(3,1) * b * DSIN(Theta_Pump) * DCOS(w) )
			END IF

		ELSE IF ( (CrystalClass .EQ. Class_4mm) .OR. ( CrystalClass .EQ. Class_6mm) ) THEN

			IF ( TypePM .EQ. 1 ) THEN

				Deff = DABS( D(3,1) * DSIN(Theta_Pump) * DCOS(-v+w) )

			ELSE

				Deff = DABS(-D(1,5) * c * DCOS(Theta_Pump) * DSIN(Phi_Pump+w) &
							+D(3,1) * a * DSIN(Theta_Pump) * DSIN(w) &
							+D(3,1) * b * DSIN(Theta_Pump) * DCOS(w) )
			END IF

		ELSE IF ( CrystalClass .EQ. Class_4_) THEN

			IF ( TypePM .EQ. 1 ) THEN

				Deff = DABS( DSIN(Theta_Pump) * ( -D(3,1) * DCOS(v+w) + D(3,6) * DSIN(v+w) ) )

			ELSE

				Deff = DABS(-D(3,6) * c * DCOS(Theta_Pump) *DCOS(Phi_Pump) * DCOS(w) &
							-D(3,1) * c * DCOS(Theta_Pump) *DCOS(Phi_Pump) * DSIN(w) &
							+D(3,1) * c * DCOS(Theta_Pump) *DSIN(Phi_Pump) * DCOS(w) &
							-D(3,6) * c * DCOS(Theta_Pump) *DSIN(Phi_Pump) * DSIN(w) &
							+D(3,1) * a * DSIN(Theta_Pump) *DSIN(w) &
							-D(3,1) * b * DSIN(Theta_Pump) *DCOS(w) &
							+D(3,6) * a * DSIN(Theta_Pump) *DCOS(w) &
							+D(3,6) * b * DSIN(Theta_Pump) *DSIN(w) ) 
			END IF

		ELSE IF ( CrystalClass .EQ. Class_6_ ) THEN

			IF ( TypePM .EQ. 1 ) THEN

				Deff = DCOS(Theta_Pump)*(D(1,1)*DCOS(Phi_Pump-v-w)-D(2,2)*DSIN(Phi_Pump-v-w))

			ELSE
				deff = -DCOS(Theta_Pump) *(DCOS(Phi_Pump)* D(1,1)* a* DSIN(w) &
						- DCOS(Phi_Pump)* D(1,1) *b* DCOS(w) - DCOS(Phi_Pump)* D(2,2)* a* DCOS(w) &
						- DCOS(Phi_Pump)* D(2,2)* b* DSIN(w) - DSIN(Phi_Pump)* D(2,2)* a* DSIN(w) &
						+ DSIN(Phi_Pump)* D(2,2)* b* DCOS(w) - DSIN(Phi_Pump)* D(1,1)* a* DCOS(w) &
						- DSIN(Phi_Pump)* D(1,1)* b* DSIN(w))
			END IF

		ELSE IF ( CrystalClass .EQ. Class_4_2m) THEN

			IF ( TypePM .EQ. 1 ) THEN

				Deff = DABS( D(3,6) * DSIN(Theta_Pump) * DSIN(v+w) )
			   
			ELSE

				Deff = DABS(- D(3,6) * ( c * DCOS(Theta_Pump) * DCOS(Phi_Pump) * DCOS(w) &
										+c * DCOS(Theta_Pump) * DSIN(Phi_Pump) * DSIN(w) &
										-a * DSIN(Theta_Pump) * DCOS(w) &
										-b * DSIN(Theta_Pump) * DSIN(w) ) )

			END IF

		ELSE IF ( CrystalClass .EQ. Class_6_2m ) THEN

			IF ( TypePM .EQ. 1 ) THEN

				Deff = DABS( D(2,2) * DCOS(Theta_Pump) * DSIN(w+v-Phi_Pump) )
				  
			ELSE

				Deff = DABS( D(2,2) * DCOS(Theta_Pump) * (DCOS(Phi_Pump) * a * DCOS(w) &
							+b * DCOS(Phi_Pump) * DSIN(w) &
							+a * DSIN(Phi_Pump) * DSIN(w) &
							-b * DSIN(Phi_Pump) * DCOS(w)) )
			   
			END IF

		END IF
	END IF

	Get_D_EFF = Deff

END FUNCTION GET_D_EFF