	PROGRAM prefdec
c**********************************************************************
c     Program decomposes clusters of overlapping reflections in a powder
c     pattern using the preferred orientation in the sample. More than
c     one powder pattern are necessary. Each with a different angle
c     between diffraction vector and vector hkl for a given reflection
c     hkl. The correction for preffered orientation according to
c     "Empirical Texture Corrections for Asymmetric Diffraction and
c      Inclined Textures", J. Appl. Cryst., Vol. 28, (1995) 247-253
c     is used.
c***********************************************************************
c     Radovan Cerny               <Radovan.Cerny@cryst.unige.ch.>
c     Laboratoire de Cristallographie
c     24, quai Ernest-Ansermet
c     CH-1211 Geneva 4, Switzerland
c     Phone: [+[41] 22] 702 64 50,  Private: [+[41] 22] 797 37 58
c     FAX  : [+[41] 22] 781 21 92,  Telex  : ch-42 11 59 siad
c************************************************************************
	PARAMETER (NDMAX=1000,MAMAX=15,NCMAX=1000,MPMAX=20,NPMAX=20)
C*****NDMAX - Max. number of pairs phiHh,y of the data
C*****MAMAX - Max. number of free variables in the texture correction
C*****NCMAX - Max. number of clusters of reflections
C*****MPMAX - Max. number of measurements of the total intensity of
C             any cluster using different angle alpha0
C*****NPMAX - Max. number of reflections in any cluster

	INTEGER hkl1(NDMAX,3),hkl2(NCMAX,NPMAX,3),ia(MAMAX),mp(NCMAX),
     *np(NCMAX),HKLtex(3),ma,ndata,j,k,l,m,nc,jump,norder

	REAL bet(NCMAX,MPMAX,NPMAX),alpha0(NCMAX,MPMAX,NPMAX),
     *b(NCMAX,MPMAX),covar(MAMAX,MAMAX),
     *mult(NCMAX,NPMAX),fwhm(NCMAX,NPMAX),yijmult(NCMAX,MPMAX,NPMAX),
     *sigb(NCMAX,MPMAX),phiHh(NPMAX),x(NPMAX),sigx(NPMAX),
     *ALPHA01(NDMAX),ALPHA02(NDMAX),sigy(NDMAX),y(NDMAX),
     *BET1(NDMAX),BET2(NDMAX),
     *latpar(6),a(MAMAX),G,sigG,wavelength,POWER,coef(15)

	CHARACTER title*80,ans,tc

	pi=3.1415926

	OPEN(UNIT=9, FILE='PATTDEC.INP')
	OPEN(UNIT=10,FILE='PREFDEC.INP')
	OPEN(UNIT=11,FILE='PREFDEC.OUT')
	OPEN(UNIT=12,FILE='SIRPOW.OUT')

	CALL prepar

	OPEN(UNIT=10,FILE='PREFDEC.INP')

	CALL input(ndata,BET1,BET2,bet,POWER,HKLtex,latpar,hkl1,hkl2,
     *y,sigy,ALPHA01,ALPHA02,ma,ia,a,nc,mp,np,alpha0,G,sigG,b,sigb,
     *title,mult,fwhm,yijmult,wavelength,NDMAX,MAMAX,NCMAX,MPMAX,NPMAX)

	CLOSE (10)

	WRITE(11,'(A80)')title
	WRITE(12,'(A80)')title

	WRITE(*,'('' Empirical (E) or spherical harmonics (S) ? '')')
	READ(*,*)tc

	WRITE(*,'('' Do you want to refine'')')
	WRITE(*,'('' the texture correction?  Y(es) or N(o)?  --> '')')
	READ(*,*)ans

	IF(ans.EQ.'Y'.OR.ans.EQ.'y') THEN
	  IF(tc.EQ.'S'.OR.tc.EQ.'s') THEN

	WRITE(*,'('' Order of spherical harmonics?  --> '')')
	READ(*,*)norder
	DO 3 i=1,norder
	 ia(i)=1.
3     CONTINUE
	WRITE(*,'('' Give JUMP --> '')')
	READ(*,*)jump
	ma=norder

	  ELSE IF(tc.EQ.'E'.OR.tc.EQ.'e') THEN

	ma=1

	  ENDIF

	CALL septex
     *(ndata,BET1,BET2,POWER,HKLtex,latpar,hkl1,y,sigy,ALPHA01,ALPHA02,
     *ma,ia,a,covar,
     *NDMAX,MAMAX,NCMAX,MPMAX,NPMAX,jump,norder,tc)

	  IF(tc.EQ.'S'.OR.tc.EQ.'s') THEN

	   DO 4 i=1,norder
	   coef(i)=a(i)
4        CONTINUE
	   write(*,*)(coef(i), i=1,norder)

	  ELSE IF(tc.EQ.'E'.OR.tc.EQ.'e') THEN

	    G=a(1)
	    sigG=SQRT(covar(1,1))

	  ENDIF

	ELSE

	  IF(tc.EQ.'S'.OR.tc.EQ.'s') THEN

	WRITE(*,'('' Order of spherical harmonics?  --> '')')
	READ(*,*)norder

	WRITE(*,'('' Give the value of '',i3,'' coef --> '')')norder
	READ(*,*)(coef(i), i=1,norder)
	WRITE(*,'('' Give JUMP --> '')')
	READ(*,*)jump

	  ELSE IF(tc.EQ.'E'.OR.tc.EQ.'e') THEN

	WRITE(*,'('' Give the value of the parameter G --> '')')
	READ(*,*)G
	WRITE(*,'('' Give the value of the st. dev. of G --> '')')
	READ(*,*)sigG
	WRITE(*,'(/,'' G= '',F9.6,4X,''sigmaG= '',F9.6,/)') G,sigG

	  ENDIF
	ENDIF

	IF(tc.EQ.'E'.OR.tc.EQ.'e') THEN
	  WRITE(11,'(/,'' G= '',F9.6,4X,''sigmaG= '',F9.6)') G,sigG
	ENDIF

	WRITE(11,'
     *(/,4X,''h'',4X,''k'',4X,''l'',12X,''Fsquare'',13X,''sigmaFsq'',
     *6X,''phiHh'')')
	WRITE(12,'
     *(3X,''h'',3X,''k'',3X,''l'',4X,''mult'',5X,''sintl'',4X,''2theta''
     *,6X,''fwhm'',3X,''Fsquare'',2X,''sigmaFsq'')')

	DO 1 i=1,nc
	CALL dectex(i,mp(i),np(i),bet,POWER,G,sigG,HKLtex,latpar,hkl2,
     *b,sigb,alpha0,phiHh,x,sigx,nc,mult,yijmult,wavelength,coef,jump,
     *NDMAX,MAMAX,NCMAX,MPMAX,NPMAX,tc)

	DO 2 j=1,np(i)
	WRITE(11,'(2X,3(I3,2X),4X,2(F16.6,2X),4X,F6.2)')
     *(hkl2(i,j,k),k=1,3),x(j),sigx(j),phiHh(j)
	dumy=rlpg(hkl2(i,j,1),hkl2(i,j,2),hkl2(i,j,3),latpar,wavelength,
     *theta)
	sintl=SIN(theta)/wavelength
	theta=theta/pi*180
	WRITE(12,'(3I4,I8,F10.5,2F10.4,F10.3,F10.4)')
     *(hkl2(i,j,k),k=1,3),mult(i,j),sintl,2*theta,fwhm(i,j),x(j),sigx(j)
2     CONTINUE
1     CONTINUE

	CLOSE (11)
	CLOSE (12)
	STOP
	END
c*********************************************************************
	SUBROUTINE prepar
c*********************************************************************
c     Subroutine uses the output from the pattern decomposition by
c     FULLPROF in the SIRPOW format (PATTDEC.INP) and prepars
c     the input file PREFDEC.INP.
C*********************************************************************
	PARAMETER (NDMAX=1000,MAMAX=15,NCMAX=1000,MPMAX=20,NPMAX=20)

	INTEGER hkl(NCMAX,3),ovl(NCMAX),nref(MPMAX),np(NCMAX),mp(NCMAX),
     *HKLtex(3),h,k,l,nreftot,nmeastot,ndata

	REAL yij(NCMAX,MPMAX),sigyij(NCMAX,MPMAX),alpha0(NCMAX,MPMAX),
     *bet(MPMAX),
     *mult(NCMAX),sintl(NCMAX),tt(NCMAX),fwhm(NCMAX),
     *latpar(6),phiHh(NCMAX),
     *incang,chi,mu,s,t,f,yy,sy,rres,dres,G,power,wavelength

	CHARACTER file(NCMAX*MPMAX)*80,geom*3,title*80,ans*1,line*80

	pi=3.1415926

c*****incang - angle (in deg) between primary beam and specimen surface,
c              has no meaning for Bragg-Brentano geometry
c*****chi    - inclination of the sample by an angle chi in Eulerian
c              cradle (in deg)
c*****bet    - angle beta of the inclination of the "inclined" texture
c              from the spinning axis. It can be different for the measured
c              powder patterns to simulate the measurements with a "precessing"
c              capillary
c*****G      - parameter of the texture correction
c*****power  - power of the texture correction
c*****geom   = BB  for Bragg-Brentano geometry
c            = SBT for Seemann-Bohlin transmission geometry
c            = SBR for Seemann-Bohlin reflection geometry
C*****ovl(i) = 0 ... reflection is not overlapped with any other reflection
C*****ovl(i) = 1 ... reflection is overlapped with the preceding reflection
C*****ovl(i) = 2 ... reflection is overlapped with an other reflection but
C                    cannot be resolved by the program PREFDEC
C*****rres - angular resolution of the pattern (2theta resolution, rec. space)
C            in multiples of FWHM
C*****dres - angular resolution of the polar vectors hkl
C*****       (orientation resolution in direct space) in deg

	READ(9,'(A80)')title
	READ(9,*)G,power,(HKLtex(i),i=1,3),(latpar(i),i=1,6),
     *         wavelength,rres,dres
	DO 1 m=1,MPMAX
	READ(9,'(A80)',END=3)line
	READ(9,'(9X,A3,3F8.0,I6,/)')geom,incang,chi,bet(m),
     *nref(m)
	write(*,*)geom,incang,chi,bet(m),nref(m)
	DO 2 i=1,nref(m)
	IF(m.EQ.1) THEN
	nreftot=nref(1)
	READ(9,'(3I4,I8,F10.5,2F10.4,2F10.0)')(hkl(i,j),j=1,3),
     *   mult(i),sintl(i),tt(i),fwhm(i),yij(i,m),sigyij(i,m)
	IF(yij(i,m).EQ.0)yij(i,m)=0.0001
C      IF(sigyij(i,m).EQ.0)sigyij(i,m)=0.1
	sigyij(i,m)=1.

	IF(geom.EQ.'BB '.OR.geom.EQ.' BB')THEN
	  alpha0(i,m)=chi
	ELSE IF(geom.EQ.'SBT')THEN
	  alpha0(i,m)=tt(i)/2+incang
	  IF(alpha0(i,m).GT.90)alpha0(i,m)=180-alpha0(i,m)
	  alpha0(i,m)=ACOS(COS(alpha0(i,m)/180*pi)*COS(chi/180*pi))*180/pi
	ELSE IF(geom.EQ.'SBR')THEN
	  alpha0(i,m)=ABS(incang-tt(i)/2)
	  alpha0(i,m)=ACOS(COS(alpha0(i,m)/180*pi)*COS(chi/180*pi))*180/pi
	ENDIF
	ELSE
	READ(9,'(3I4,I8,F10.5,2F10.4,2F10.0)')h,k,l,mu,s,t,f,yy,sy
	IF(yy.EQ.0)yy=0.0001
c      IF(sy.EQ.0)sy=0.1
	sy=1.

	IF(geom.EQ.'BB '.OR.geom.EQ.' BB')THEN
	  al0=chi
	ELSE IF(geom.EQ.'SBT')THEN
	  al0=t/2+incang
	  IF(al0.GT.90)al0=180-al0
	  al0=ACOS(COS(al0/180*pi)*COS(chi/180*pi))*180/pi
	ELSE IF(geom.EQ.'SBR')THEN
	  al0=ABS(incang-t/2)
	  al0=ACOS(COS(al0/180*pi)*COS(chi/180*pi))*180/pi
	ENDIF
	DO 4 j=1,nreftot
	IF(h.EQ.hkl(j,1).AND.k.EQ.hkl(j,2).AND.l.EQ.hkl(j,3))THEN
	  yij(j,m)=yy
	  sigyij(j,m)=sy
	  alpha0(j,m)=al0
	  GOTO 5
	ENDIF
4     CONTINUE
	  nreftot=nreftot+1
	  hkl(nreftot,1)=h
	  hkl(nreftot,2)=k
	  hkl(nreftot,3)=l
	  mult(nreftot)=mu
	  sintl(nreftot)=s
	  tt(nreftot)=t
	  fwhm(nreftot)=f
	  yij(nreftot,m)=yy
	  sigyij(nreftot,m)=sy
	  alpha0(nreftot,m)=al0
5     CONTINUE
	ENDIF
2     CONTINUE

1     CONTINUE
3     nmeastot=m-1
	ovl(1)=0

	DO 6 i=2,nreftot
	ovl(i)=0
	IF(ABS(tt(i)-tt(i-1)).LE.(rres*fwhm(i-1))) ovl(i)=1
6     CONTINUE

	DO 19 i=1,nreftot
	IF(ovl(i).EQ.0)THEN
	  DO 20 j=i,nreftot
	  phiHh(j)=angPhiHh(hkl(j,1),hkl(j,2),hkl(j,3),
     *  HKLtex(1),HKLtex(2),HKLtex(3),
     *  latpar(1),latpar(2),latpar(3),latpar(4),latpar(5),latpar(6))
	  IF(ovl(j+1).EQ.0)THEN
	  ii=j
	  GOTO 22
	  ENDIF
20      CONTINUE
	  ii=nreftot
22      DO 21 k=i,ii-1
	     DO 21 l=i+1,ii
	     IF((ABS(phiHh(k)-phiHh(l)).LT.dres).AND.(k.NE.l)
     *     .AND.(ABS(tt(k)-tt(l)).LT.(l-k)*rres*fwhm(k)))THEN
	     ovl(k)=2
c           ovl(l)=2
	     ENDIF
21      CONTINUE
	ENDIF
19    CONTINUE

	DO 23 i=1,nreftot-1
	IF(ovl(i).EQ.2)THEN
	  DO 24 j=i+1,nreftot
	    IF(ovl(j).NE.2)THEN
	    ovl(j)=0
	    GOTO 23
	    ENDIF
24      CONTINUE
	ENDIF
23    CONTINUE

	ovl(nreftot+1)=0
	ndata=0
	DO 10 i=1,nreftot
	IF(ovl(i).EQ.0.AND.(ovl(i+1).EQ.0.OR.ovl(i+1).EQ.2))THEN
	  DO 8 j=1,nmeastot-1
	     DO 8 k=j+1,nmeastot
	     IF(yij(i,j).NE.0.AND.yij(i,k).NE.0.AND.j.NE.k)THEN
		  ndata=ndata+1
		  WRITE(file(ndata),'(3(I3,1X),2(F12.6,1X),4(F8.4,1X))')
     *        (hkl(i,l),l=1,3),yij(i,j)/yij(i,k),
     *        SQRT(1/(yij(i,k)**2)*sigyij(i,j)**2+
     *             yij(i,j)**2/yij(i,k)**4*sigyij(i,k)**2),
     *        alpha0(i,j),alpha0(i,k),bet(j),bet(k)
	     ENDIF
8          CONTINUE
	ENDIF
	IF(ndata.GE.40)GOTO 26
10    CONTINUE

26    WRITE(10,'(A80)')title
	WRITE(10,'(I3,1X,F4.1,1X,3(I3,1X),3(F7.4,1X),3(F6.2,1X),F6.4)')
     *ndata,power,(HKLtex(i),i=1,3),(latpar(i),i=1,6),wavelength

	DO 9 i=1,ndata
	READ(file(i),'(A80)')line
	WRITE(10,'(A80)')line
9     CONTINUE
	WRITE(10,'(''1 '',/,''1 '',/,F8.4)')G

	nc=0
	n=1
	DO 11 i=1,nreftot
	IF(ovl(i).EQ.1)GOTO 11
	   nsave=n
	   nc=nc+1
25       mp(nc)=0
	   n=nsave
	   DO 14 k=1,nmeastot
	   IF(yij(i,k).NE.0)THEN
	   b=yij(i,k)*mult(i)
	   sigb=sigyij(i,k)*mult(i)
	   DO 12 j=i+1,nreftot
	   IF(ovl(j).EQ.1)THEN
	     b=b+yij(j,k)*mult(j)
	     sigb=SQRT(sigb**2+(sigyij(j,k)*mult(j))**2)
	   ELSE
	     GOTO 13
	   ENDIF
12       CONTINUE
13       np(nc)=j-1-i+1
	   mp(nc)=mp(nc)+1
	   WRITE(file(n),*)b,sigb
	   n=n+1
		DO 15 l=i,j-1
	   WRITE(file(n),'(3(I3,1X),2(F6.2,1X),F3.0,1X,F7.5,1X,2(F9.4,1X)
     *      ,1X,F12.4,1X,I2)')
     *      (hkl(l,m),m=1,3),alpha0(l,k),bet(k),mult(l),fwhm(l),tt(l),
     *      phiHh(l),yij(l,k),ovl(l)
		n=n+1
15          CONTINUE
	   ENDIF
14       CONTINUE

	   IF(np(nc).GT.mp(nc))THEN
	     ovl(INT((j-1+i+1)/2))=0
	     GOTO 25
	   ENDIF

11    CONTINUE
	n=1
	WRITE(10,'(I4)')nc
	DO 16 i=1,nc
	   WRITE(10,'(2(I4,1X))')mp(i),np(i)
	   DO 17 j=1,mp(i)
	   READ(file(n),'(A80)')line
	   n=n+1
	   WRITE(10,'(A80)')line
		 DO 18 k=1,np(i)
		    READ(file(n),'(A80)')line
		    n=n+1
		    WRITE(10,'(A80)')line
18           CONTINUE
17       CONTINUE
16    CONTINUE

	CLOSE(10)
	PAUSE
	RETURN
	END
c*********************************************************************
	SUBROUTINE input(ndata,BET1,BET2,bet,POWER,HKLtex,latpar,hkl1,
     *hkl2,y,sigy,ALPHA01,ALPHA02,ma,ia,a,nc,mp,np,alpha0,G,sigG,b,
     *sigb,title,mult,fwhm,yijmult,wavelength,
     *NDMAX,MAMAX,NCMAX,MPMAX,NPMAX)
c*********************************************************************
C     Subroutine reads the input file PREFDEC.INP
c**********************************************************************
	INTEGER ma,ndata,j,k,l,m,ia(MAMAX),nc,mp(NCMAX),np(NCMAX),
     *HKLtex(3),hkl1(NDMAX,3),hkl2(NCMAX,NPMAX,3)

	REAL latpar(6),a(MAMAX),sigy(NDMAX),y(NDMAX),ALPHA01(NDMAX),
     *ALPHA02(NDMAX),BET1(NDMAX),BET2(NDMAX),bet(NCMAX,MPMAX,NPMAX),
     *POWER,alpha0(NCMAX,MPMAX,NPMAX),G,sigG,
     *b(NCMAX,MPMAX),sigb(NCMAX,MPMAX),mult(NCMAX,NPMAX),
     *fwhm(NCMAX,NPMAX),yijmult(NCMAX,MPMAX,NPMAX)
	CHARACTER title*80

C*****ndata  - number of pairs phiHh,y of the data
c*****BET    - angle beta of the inclination of the "inclined" texture
c              from the spinning axis. It can be different for the measured
c              powder patterns to simulate the measurements with a "precessing"
c              capillary
C*****POWER  - power of the exponent in the texture correction
C*****HKLtex - hkls of preferrentially oriented plane
C*****G      - texture parameter
C*****hkl1   - hkls of planes used for the refinement of texture correction
	READ(10,'(A80)')title
	WRITE(*,*)title
	READ(10,*)ndata,POWER,(HKLtex(i),i=1,3),(latpar(i),i=1,6),
     *          wavelength

	DO 1 i=1,ndata

C*****y(i), sigy(i) - ratio between the intensities of reflection hkl
C*****               measured at the angle ALPHA01 and ALPHA02
C*****ALPHA01 and ALPHA02 - angles between the diffraction vector and
C*****                      the spinning axis
	READ(10,*)(hkl1(i,j),j=1,3),y(i),sigy(i),ALPHA01(i),ALPHA02(i),
     *          BET1(i),BET2(i)
1     CONTINUE

C*****ma - number of free variables in the texture correction
	READ(10,*)ma

C*****ia(i) = 0 - i's parameter is fixed
C*****ia(i) = 1 - i's parameter is refined
	READ(10,*)(ia(i),i=1,ma)

C*****a(i) - starting estimate for the i's parameter
	READ(10,*)(a(i),i=1,ma)

C*****nc    - number of clusters of reflections
C*****mp(i) - number of measurements of the total intensity of the
C             i-th cluster using different angle alpha0
C*****np(i) - number of reflections in the i-th cluster


	READ(10,*)nc

	DO 4 k=1,nc
	READ(10,*)mp(k),np(k)

	DO 3 i=1,mp(k)
C*****b(k,i) - total intensity of k-th cluster measured with i-th angle alpha0
	READ(10,*)b(k,i),sigb(k,i)

	DO 2 j=1,np(k)
C*****hkl2 - hkls of reflections to be separated
C*****mult - multiplicity of the reflection
C*****alpha0(k,i,j) - angle between the diffraction vector and
C                 the spinning axis
C*****yijmult - Fsqr of j-th reflection in k-th cluster multiplied by the
C               reflection multipplicity
	READ(10,*)(hkl2(k,j,l),l=1,3),alpha0(k,i,j),bet(k,i,j),
     *          mult(k,j),fwhm(k,j),dummy,dummy,yijmult(k,i,j)
2     CONTINUE
3     CONTINUE
4     CONTINUE

	RETURN
	END
c*******************************************************************
	SUBROUTINE septex
     *(ndata,BET1,BET2,POWER,HKLtex,latpar,hkl1,y,sigy,ALPHA01,ALPHA02,
     *ma,ia,a,covar,
     *NDMAX,MAMAX,NCMAX,MPMAX,NPMAX,jump,norder,tc)
c*******************************************************************
c     Subroutine refines the parameters of the texture correction from
c     the ratios of intensities of the reflections hkl measured at
c     different angles ALPHA0.

c     Radovan Cerny               <Radovan.Cerny@cryst.unige.ch.>
c     Laboratoire de Cristallographie    
c     24, quai Ernest-Ansermet              
c     CH-1211 Geneva 4, Switzerland                            
c     Phone: [+[41] 22] 702 64 50,  Private: [+[41] 22] 797 37 58
c     FAX  : [+[41] 22] 781 21 92,  Telex  : ch-42 11 59 siad
c*******************************************************************
	INTEGER ma,nca,ndata,MMAX,HKLtex(3),hkl1(NDMAX,3),jump,norder
	REAL latpar(6),alamda,chisq,a(MAMAX),astart(:),asave(:),
     *alpha(:,:),covar(MAMAX,MAMAX),sigy(NDMAX),phiHh(:),y(NDMAX),
     *ALPHA01(NDMAX),ALPHA02(NDMAX),BET1(NDMAX),BET2(NDMAX),POWER
	PARAMETER (MMAX=15)
CU    USES mrqmin
	INTEGER j,k,l,m,mfit,ia(MAMAX)
	ALLOCATABLE astart,asave,alpha,phiHh
	REAL ochisq,atry(MMAX),beta(MMAX),da(MMAX)
	CHARACTER tc

	nca=0

	DO 2 i=1,ma
	IF(ia(i).EQ.1)nca=nca+1
2     CONTINUE
	ALLOCATE (astart(ma),asave(ma),alpha(nca,nca),phiHh(ndata))

	DO 3 i=1,ndata
	phiHh(i)=angPhiHh(hkl1(i,1),hkl1(i,2),hkl1(i,3),
     *HKLtex(1),HKLtex(2),HKLtex(3),
     *latpar(1),latpar(2),latpar(3),latpar(4),latpar(5),latpar(6))
3     CONTINUE

	DO 10 k=1,ma
	astart(k)=a(k)
10    CONTINUE
	l=1
	alamda=-1.

9     DO 4 i=l,l+1
	DO 6 k=1,ma
	asave(k)=a(k)
6     CONTINUE

	WRITE(11,'
     *(''  PHIHH     ALPHA01  DERIV1   ALPHA02  '',
     *'' DERIV2   YFIT    DYDA     Y(i)'')')

	CALL mrqmin(phiHh,y,sigy,ndata,a,ia,ma,covar,alpha,nca,chisq,
     *alamda,ALPHA01,ALPHA02,BET1,BET2,POWER,
     *NDMAX,MAMAX,NCMAX,MPMAX,NPMAX,latpar,hkl1,jump,norder,tc)

	WRITE(6,'(/'' Cycle no. '',I3)')i
	WRITE(6,'(E10.4,2X,E8.4,10X,15E8.4)')
     *           alamda,chisq,(a(j),j=1,ma)
	WRITE(11,'(/'' Cycle no. '',I3)')i
	WRITE(11,'(E10.4,2X,E8.4,10X,15E8.4)')
     *           alamda,chisq,(a(j),j=1,ma)
4     CONTINUE
	DO 8 i=1,ma
	IF((ABS(asave(i)-a(i)).GT.0.0001.OR.a(i).EQ.astart(i))
     *   .AND.l.LE.8) THEN
	l=l+2
	GOTO 9
	END IF
8     CONTINUE
	alamda=0.

	CALL mrqmin(phiHh,y,sigy,ndata,a,ia,ma,covar,alpha,nca,chisq,
     *alamda,ALPHA01,ALPHA02,BET1,BET2,POWER,
     *NDMAX,MAMAX,NCMAX,MPMAX,NPMAX,latpar,hkl1,jump,norder,tc)

	WRITE(6,'(/''  Lambda    Chisqr                 Parameters'')')
	WRITE(6,'(E10.4,1X,E10.4,8X,15E10.4)')
     *           alamda,chisq,(a(j),j=1,ma)
	WRITE(6,'(''                                      E.S.D.'')')
	WRITE(6,'(30X,15E10.4)')(SQRT(ABS(covar(k,k))),k=1,ma)
	WRITE(11,'(/''  Lambda    Chisqr                 Parameters'')')
	WRITE(11,'(E10.4,1X,E10.4,8X,15E10.4)')
     *           alamda,chisq,(a(j),j=1,ma)
	WRITE(11,'(''                                      E.S.D.'')')
	WRITE(11,'(30X,15E10.4)')(SQRT(ABS(covar(k,k))),k=1,ma)
	PAUSE
	DEALLOCATE (astart,asave,alpha,phiHh)
	RETURN
	END
c*******************************************************************
	SUBROUTINE mrqmin(x,y,sigy,ndata,a,ia,ma,covar,alpha,nca,chisq,
     *alamda,ALPHA01,ALPHA02,BET1,BET2,POWER,
     *NDMAX,MAMAX,NCMAX,MPMAX,NPMAX,latpar,hkl1,jump,norder,tc)
c*******************************************************************
	INTEGER ma,nca,ndata,ia(MAMAX),MMAX,hkl1(NDMAX,3),jump,norder
	REAL alamda,chisq,a(MAMAX),alpha(nca,nca),covar(MAMAX,MAMAX),
     *sigy(NDMAX),x(ndata),y(NDMAX),ALPHA01(NDMAX),ALPHA02(NDMAX),
     *BET1(NDMAX),BET2(NDMAX),latpar(6),
	PARAMETER (MMAX=15)
CU    USES covsrt,gaussj,mrqcof
	INTEGER j,k,l,m,mfit
	REAL ochisq,atry(MMAX),beta(MMAX),da(MMAX)
	CHARACTER tc

	SAVE ochisq,atry,beta,da,mfit

	if(alamda.lt.0.)then
        mfit=0
        do 11 j=1,ma
          if (ia(j).ne.0) mfit=mfit+1
11      continue
        alamda=0.001
	  call mrqcof(x,y,sigy,ndata,a,ia,ma,alpha,beta,nca,chisq,
     *  ALPHA01,ALPHA02,BET1,BET2,POWER,
     *  NDMAX,MAMAX,NCMAX,MPMAX,NPMAX,latpar,hkl1,jump,norder,tc)

        ochisq=chisq
        do 12 j=1,ma
          atry(j)=a(j)
12      continue
      endif
      j=0
      do 14 l=1,ma
        if(ia(l).ne.0) then
          j=j+1
          k=0
          do 13 m=1,ma
            if(ia(m).ne.0) then
              k=k+1
              covar(j,k)=alpha(j,k)
            endif
13        continue
          covar(j,j)=alpha(j,j)*(1.+alamda)
          da(j)=beta(j)
        endif
14    continue

	call gaussj(covar,mfit,nca,da,1,1,
     *NDMAX,MAMAX,NCMAX,MPMAX,NPMAX)

      if(alamda.eq.0.)then
	  call covsrt(covar,nca,ma,ia,mfit,
     *NDMAX,MAMAX,NCMAX,MPMAX,NPMAX)

        return
      endif
      j=0
      do 15 l=1,ma
        if(ia(l).ne.0) then
          j=j+1
          atry(l)=a(l)+da(j)
        endif
15    continue
	call mrqcof(x,y,sigy,ndata,atry,ia,ma,covar,da,nca,chisq,
     *ALPHA01,ALPHA02,BET1,BET2,POWER,
     *NDMAX,MAMAX,NCMAX,MPMAX,NPMAX,latpar,hkl1,jump,norder,tc)

	if(chisq.lt.ochisq)then
        alamda=0.1*alamda
        ochisq=chisq
        j=0
        do 17 l=1,ma
          if(ia(l).ne.0) then
            j=j+1
            k=0
            do 16 m=1,ma
              if(ia(m).ne.0) then
                k=k+1
                alpha(j,k)=covar(j,k)
              endif
16          continue
            beta(j)=da(j)
            a(l)=atry(l)
          endif
17      continue
	else
	  alamda=10.*alamda
	  chisq=ochisq
	endif
      return
	END
c***********************************************************************
	SUBROUTINE mrqcof(x,y,sigy,ndata,a,ia,ma,alpha,beta,nalp,chisq,
     *ALPHA01,ALPHA02,BET1,BET2,POWER,
     *NDMAX,MAMAX,NCMAX,MPMAX,NPMAX,latpar,hkl1,jump,norder,tc)
c***********************************************************************
	INTEGER ma,nalp,ndata,ia(MAMAX),MMAX,hkl3(3),hkl1(NDMAX,3),
     *        jump,norder
	REAL chisq,a(MAMAX),alpha(nalp,nalp),beta(ma),sigy(NDMAX),
     *x(ndata),y(NDMAX),ALPHA01(NDMAX),ALPHA02(NDMAX),latpar(6)
     *BET1(NDMAX),BET2(NDMAX)
	PARAMETER (MMAX=15)
      INTEGER mfit,i,j,k,l,m
	REAL dy,sig2i,wt,ymod,dyda(MMAX)
	CHARACTER tc

      mfit=0
      do 11 j=1,ma
        if (ia(j).ne.0) mfit=mfit+1
11    continue
      do 13 j=1,mfit
        do 12 k=1,j
          alpha(j,k)=0.
12      continue
        beta(j)=0.
13    continue
      chisq=0.
	do 16 i=1,ndata
	 hkl3(1)=hkl1(i,1)
	 hkl3(2)=hkl1(i,2)
	 hkl3(3)=hkl1(i,3)
	 call funcs(x(i),a,ymod,dyda,ma,ALPHA01(i),ALPHA02(i),
     * BET1(i),BET2(i),POWER,
     * NDMAX,MAMAX,NCMAX,MPMAX,NPMAX,y(i),latpar,hkl3,jump,norder,tc)
	  sig2i=1./(sigy(i)*sigy(i))
        dy=y(i)-ymod
        j=0
        do 15 l=1,ma
          if(ia(l).ne.0) then
            j=j+1
            wt=dyda(l)*sig2i
            k=0
            do 14 m=1,l
              if(ia(m).ne.0) then
                k=k+1
                alpha(j,k)=alpha(j,k)+wt*dyda(m)
              endif
14          continue
            beta(j)=beta(j)+dy*wt
          endif
15      continue
        chisq=chisq+dy*dy*sig2i
16    continue
      do 18 j=2,mfit
        do 17 k=1,j-1
          alpha(k,j)=alpha(j,k)
17      continue
18    continue
      return
	END
c***********************************************************************
	SUBROUTINE covsrt(covar,npc,ma,ia,mfit,
     *NDMAX,MAMAX,NCMAX,MPMAX,NPMAX)
c***********************************************************************
	INTEGER ma,mfit,npc,ia(MAMAX)
	REAL covar(MAMAX,MAMAX)
      INTEGER i,j,k
      REAL swap
      do 12 i=mfit+1,ma
        do 11 j=1,i
          covar(i,j)=0.
          covar(j,i)=0.
11      continue
12    continue
      k=mfit
      do 15 j=ma,1,-1
        if(ia(j).ne.0)then
          do 13 i=1,ma
            swap=covar(i,k)
            covar(i,k)=covar(i,j)
            covar(i,j)=swap
13        continue
          do 14 i=1,ma
            swap=covar(k,i)
            covar(k,i)=covar(j,i)
            covar(j,i)=swap
14        continue
          k=k-1
        endif
15    continue
      return
	END
c************************************************************************
	SUBROUTINE gaussj(a,n,np,b,m,mp,
     *NDMAX,MAMAX,NCMAX,MPMAX,NPMAX)
c************************************************************************
      INTEGER m,mp,n,np,NMAX
	REAL a(MAMAX,MAMAX),b(np,mp)
      PARAMETER (NMAX=50)
      INTEGER i,icol,irow,j,k,l,ll,indxc(NMAX),indxr(NMAX),ipiv(NMAX)
      REAL big,dum,pivinv
      do 11 j=1,n
        ipiv(j)=0
11    continue
      do 22 i=1,n
        big=0.
        do 13 j=1,n
          if(ipiv(j).ne.1)then
            do 12 k=1,n
              if (ipiv(k).eq.0) then
                if (abs(a(j,k)).ge.big)then
                  big=abs(a(j,k))
                  irow=j
                  icol=k
                endif
              else if (ipiv(k).gt.1) then
                pause 'singular matrix in gaussj'
              endif
12          continue
          endif
13      continue
        ipiv(icol)=ipiv(icol)+1
        if (irow.ne.icol) then
          do 14 l=1,n
            dum=a(irow,l)
            a(irow,l)=a(icol,l)
            a(icol,l)=dum
14        continue
          do 15 l=1,m
            dum=b(irow,l)
            b(irow,l)=b(icol,l)
            b(icol,l)=dum
15        continue
        endif
        indxr(i)=irow
        indxc(i)=icol
        if (a(icol,icol).eq.0.) pause 'singular matrix in gaussj'
        pivinv=1./a(icol,icol)
        a(icol,icol)=1.
        do 16 l=1,n
          a(icol,l)=a(icol,l)*pivinv
16      continue
        do 17 l=1,m
          b(icol,l)=b(icol,l)*pivinv
17      continue
        do 21 ll=1,n
          if(ll.ne.icol)then
            dum=a(ll,icol)
            a(ll,icol)=0.
            do 18 l=1,n
              a(ll,l)=a(ll,l)-a(icol,l)*dum
18          continue
            do 19 l=1,m
              b(ll,l)=b(ll,l)-b(icol,l)*dum
19          continue
          endif
21      continue
22    continue
      do 24 l=n,1,-1
        if(indxr(l).ne.indxc(l))then
          do 23 k=1,n
            dum=a(k,indxr(l))
            a(k,indxr(l))=a(k,indxc(l))
            a(k,indxc(l))=dum
23        continue
        endif
24    continue
      return
	END
c****************************************************************************
C  (C) Copr. 1986-92 Numerical Recipes Software Y"1i41)jD#ks<,1'i0S%.
c****************************************************************************
	SUBROUTINE funcs(PHIHH,a,yfit,dyda,ma,ALPH01,ALPH02,BET1,BET2,
     *                 POWER,NDMAX,MAMAX,NCMAX,MPMAX,NPMAX,yobs,
     *                 latpar,hkl3,jump,norder,tc)
c*******************************************************************
	INTEGER ma,MMAX,hkl3(3),jump,norder
	PARAMETER (MMAX=15)
	REAL a(MAMAX),yfit,dyda(MMAX),PHIHH,ALPH01,ALPH02,BET1,BET2,
     *     POWER,yobs,latpar(6),texder(MMAX)
	CHARACTER tc

	DEG2RAD=3.1415927/180

C*****All angles ARE in DEG.**

	IF(tc.EQ.'E'.OR.tc.EQ.'e') THEN

C*****Using NEW TEXTURE CORRECTION FOR INCLINED TEXTURE AND ASYMMETRIC DIFFRACTION
	P1=TEXTUR(PHIHH*DEG2RAD,ALPH01*DEG2RAD,BET1,POWER,a(1),1)
	P2=TEXTUR(PHIHH*DEG2RAD,ALPH02*DEG2RAD,BET2,POWER,a(1),1)
	write(11,*)phihh,alph01,alph02,bet,power,a(1),p1,p2
	yfit=P1/P2
	d1=TEXTUR(PHIHH*DEG2RAD,ALPH01*DEG2RAD,BET1,POWER,a(1),2)
	d2=TEXTUR(PHIHH*DEG2RAD,ALPH02*DEG2RAD,BET2,POWER,a(1),2)
	dyda(1)=d1/P2-P1/P2**2*d2
	WRITE(11,'(10F8.3)') PHIHH,ALPH01,BET1,d1,ALPH02,BET2,d2,
     *                    yfit,dyda(1),yobs
C*****End of using NEW TEXTURE CORRECTION FOR INCLINED TEXTURE AND ASYMMETRIC DIFFRACTION

	ELSE IF(tc.EQ.'S'.OR.tc.EQ.'s') THEN

C*****Using texture correction with symmtrized spherical harmonics.
	geom=COS(ALPH01*DEG2RAD)
	CALL texsphh(latpar,hkl3(1),hkl3(2),hkl3(3),
     *  jump,geom,P1,a,texder,norder)
	do 2 i=1,15
	dyda(i)=texder(i)
	dumy=texder(1)
2     continue
	geom=COS(ALPH02*DEG2RAD)
	CALL texsphh(latpar,hkl3(1),hkl3(2),hkl3(3),
     *  jump,geom,P2,a,texder,norder)
	yfit=P1/P2
	do 1 i=1,15
	dyda(i)=dyda(i)/P2-P1/P2**2*texder(i)
1     continue
	WRITE(11,'(10F8.3)') PHIHH,ALPH01,BET1,DUMY,ALPH02,BET2,TEXDER(2),
     *                    yfit,dyda(1),yobs
C*****Using texture correction with symmtrized spherical harmonics.

	ENDIF

	RETURN
	END
c***********************************************************************
	SUBROUTINE dectex
     *(ic,mp,np,bet,POWER,G,sigG,HKLtex,latpar,hkl2,b,sigb,alpha0,phiHh,
     *x,sigx,nc,mult,yijmult,wavelength,coef,jump,
     *NDMAX,MAMAX,NCMAX,MPMAX,NPMAX,tc)
c***********************************************************************
c     Subroutine finds the vaules of F2 of N overlapping reflections
c     in a cluster from M measurements of the total intensity of
c     this cluster with different angles alpha0. It needs to know
c     already the texture function PHKL (FUNCTION TEXTUR) found
c     by the subroutine septex.

c     Radovan Cerny               <Radovan.Cerny@cryst.unige.ch.>
c     Laboratoire de Cristallographie    
c     24, quai Ernest-Ansermet              
c     CH-1211 Geneva 4, Switzerland                            
c     Phone: [+[41] 22] 702 64 50,  Private: [+[41] 22] 797 37 58
c     FAX  : [+[41] 22] 781 21 92,  Telex  : ch-42 11 59 siad
c*******************************************************************
	INTEGER m,mp,n,np,nc,scr,kbd,HKLtex(3),
     *hkl2(NCMAX,NPMAX,3),jump,norder

	REAL latpar(6),a(:,:),siga(:,:),u(:,:),w(:),v(:,:),
     *b(NCMAX,MPMAX),sigb(NCMAX,MPMAX),bb(:),x(NPMAX),sigx(NPMAX),
     *phiHh(NPMAX),alpha0(NCMAX,MPMAX,NPMAX),bet(NCMAX,MPMAX,NPMAX),
     *POWER,G,sigG,bcalc,mult(NCMAX,NPMAX),yijmult(NCMAX,MPMAX,NPMAX),
     *wavelength,coef(15),texturecor,texder(15)
	CHARACTER tc

	PARAMETER (NMAX=1000)

	ALLOCATABLE a,siga,u,w,v,bb
CU    USES svdcmp,svbksp,textur
	pi=3.1415927
	ALLOCATE (a(mp,np),siga(mp,np),u(mp+np,np),w(np),v(np,np),
     *         bb(mp+np))

C*****ic       - serial number of a cluster
C*****phiHh(j) - angle phiH,h (in DEG) between the plane hkl and
C                the textured plane HKL

	DO 3 i=1,mp
	bb(i)=b(ic,i)
	  DO 4 j=1,np
	  phiHh(j)=angPhiHh(hkl2(ic,j,1),hkl2(ic,j,2),hkl2(ic,j,3),
     *  HKLtex(1),HKLtex(2),HKLtex(3),
     *  latpar(1),latpar(2),latpar(3),latpar(4),latpar(5),latpar(6))

	IF(tc.EQ.'S'.OR.tc.EQ.'s') THEN

C*****Using texture correction with symmetrized spherical harmonics.
	geom=COS(alpha0(ic,i,j)/180*pi)
	CALL texsphh(latpar,hkl2(ic,j,1),hkl2(ic,j,2),hkl2(ic,j,3),
     *  jump,geom,texturecor,coef,texder,norder)
	  a(i,j)=texturecor*mult(ic,j)
C*****End of using texture correction with symmetrized spherical harmonics.

	ELSE IF(tc.EQ.'E'.OR.tc.EQ.'e') THEN

C*****Using NEW TEXTURE CORRECTION FOR INCLINED TEXTURE AND ASYMMETRIC DIFFRACTION
	  a(i,j)=
     *  TEXTUR(phiHh(j)/180*pi,alpha0(ic,i,j)/180*pi,bet(ic,i,j),
     *         POWER,G,1)*mult(ic,j)
C*****End of using NEW TEXTURE CORRECTION FOR INCLINED TEXTURE AND ASYMMETRIC DIFFRACTION

	ENDIF

c     *  *rlpg(hkl2(ic,j,1),hkl2(ic,j,2),hkl2(ic,j,3),latpar,wavelength,
c     *   theta)
C*****phiHh/180*pi and alpha0/180*pi in RAD
C*****BET in DEG
	   u(i,j)=a(i,j)
4        CONTINUE
3     CONTINUE

	DO 20 j=1,np
	 sumu=0.
	 sumuyijm=0.
	 DO 21 i=1,mp
	  sumu=sumu+u(i,j)
	  sumuyijm=sumuyijm+u(i,j)*yijmult(ic,i,j)
21    CONTINUE
	 u(mp+j,j)=sumu/sumuyijm
	 bb(mp+j)=1.
20    CONTINUE

C  Hard link to the Le Bail method :
C      GOTO 15

C  This is the single value decomposition method :
	CALL svdcmp(u,mp+np,np,mp+np,np,w,v)
	wmax=0.
	DO 5 j=1,np
	   IF(w(j).GT.wmax)wmax=w(j)
5     CONTINUE
	wmin=wmax*1.0E-06
	DO 6 j=1,np
	IF(w(j).LT.wmin)w(j)=0.
6     CONTINUE
	CALL svbksb(u,w,v,mp+np,np,mp+np,np,bb,x,
     *NDMAX,MAMAX,NCMAX,MPMAX,NPMAX)

c      DO 10 i=1,mp
c      bcalc=0.
c      DO 13 j=1,np
c         bcalc=bcalc+a(i,j)*x(j)
c13    CONTINUE
c      WRITE(6,'('' No.  calc. right side        obs. right side'')')
c      WRITE(6,'(1X,I3,5X,F12.6,12X,F12.6)')i,bcalc,b(ic,i)
c10    CONTINUE

	DO 8 i=1,mp
	   DO 9 j=1,np
	   siga(i,j)=0.
9        CONTINUE
8     CONTINUE

	DO 7 k=1,np
	sigxmax=1.E+22
	sigx(k)=sigxmax
	DO 11 i=1,mp
	sigxk=0.
	   DO 12 j=1,np
	   sigxk=sigxk+x(j)*siga(i,j)
12        CONTINUE
	sigxk=sigxk+sigb(ic,i)
	sigxk=sigxk/a(i,k)
	IF((sigxk.LT.sigx(k)).AND.(sigxk.NE.0.))sigx(k)=sigxk
11    CONTINUE

c      WRITE(6,'(''      reflection         Fsquare'')')
c      WRITE(6,'(6X,I3,15X,F12.6,1X,''+-'',F12.6)')k,x(k),sigx(k)
7     CONTINUE
C  This is the end  of the single value decomposition method :

C  This is the iteration relaxation method :
C14    CALL iterelax(mp,np,u,bb,x,NPMAX)
C  This is the end of the iteration relaxation method :

C  This is Le Bail method :
C15    CALL lebail(mp,np,u,bb,x,NPMAX)
C  This is the end of Le Bail method :

	DEALLOCATE (a,siga,u,w,v,bb)
	RETURN
	END
c***********************************************************************
	SUBROUTINE lebail(m,n,u,bb,x,NPMAX)
c***********************************************************************
c     Le Bail method.

c     Radovan Cerny               <Radovan.Cerny@cryst.unige.ch.>
c     Laboratoire de Cristallographie
c     24, quai Ernest-Ansermet
c     CH-1211 Geneva 4, Switzerland
c     Phone: [+[41] 22] 702 64 50,  Private: [+[41] 22] 797 37 58
c     FAX  : [+[41] 22] 781 21 92,  Telex  : ch-42 11 59 siad
c*******************************************************************
	INTEGER m,n
	REAL u(m,n),bb(m),sumu,sumuxj,xoldsum,xnewsum,sumuxi,
     *     x(NPMAX),xold(:),xnew(:),xnewj(:,:),xoldj(:,:)

	ALLOCATABLE xold,xnew,xoldj,xnewj
	ALLOCATE (xold(NPMAX),xnew(NPMAX),xoldj(m,n),xnewj(m,n))


	do 1 j=1,m
	 do 1 i=1,n
	 xoldj(j,i)=bb(j)/n/u(j,i)
1     continue

	do 3 i=1,n
	sumu=0.
	sumuxj=0.
	 do 4 j=1,m
	  sumu=sumu+u(j,i)
	  sumuxj=sumuxj+u(j,i)*xoldj(j,i)
4     continue
	 xold(i)=sumuxj/sumu
3     continue

8     do 5 j=1,m
	 do 5 i=1,n
	sumuxi=0.
	 do 6 ii=1,n
	    sumuxi=sumuxi+u(j,ii)*xold(ii)
6      continue
	 xnewj(j,i)=xoldj(j,i)*(bb(j)/sumuxi)
c       if(xnewj(j,i).LT.0)xnewj(j,i)=0.
5     continue

	xoldsum=0.
	xnewsum=0.
	do 13 i=1,n
	sumu=0.
	sumuxj=0.
	 do 14 j=1,m
	  sumu=sumu+u(j,i)
	  sumuxj=sumuxj+u(j,i)*xnewj(j,i)
14    continue
	 xnew(i)=sumuxj/sumu
	 xoldsum=xoldsum+xold(i)
	 xnewsum=xnewsum+xnew(i)
13    continue

c      write(*,*)(xold(i), i=1,n)
c      write(*,*)(xnew(i), i=1,n)
c      pause

	if(ABS((xoldsum-xnewsum)/xnewsum).GT.0.05)then
	 do 15 j=1,m
	  do 15 i=1,n
	  xoldj(j,i)=xnewj(j,i)
	  xold(i)=xnew(i)
15     continue
	 goto 8
	endif

	do 9 i=1,n
	 x(i)=xnew(i)
9     continue

	DEALLOCATE (xold,xnew,xoldj,xnewj)
	RETURN
	END
c***********************************************************************
	SUBROUTINE iterelax(m,n,u,bb,x,NPMAX)
c***********************************************************************
c     Solve a system of m equations for n unknowns (m>=n) by the relaxation
c     method (Rektorys page 951).

c     Radovan Cerny               <Radovan.Cerny@cryst.unige.ch.>
c     Laboratoire de Cristallographie
c     24, quai Ernest-Ansermet
c     CH-1211 Geneva 4, Switzerland
c     Phone: [+[41] 22] 702 64 50,  Private: [+[41] 22] 797 37 58
c     FAX  : [+[41] 22] 781 21 92,  Telex  : ch-42 11 59 siad
c*******************************************************************
	INTEGER m,n,p,pp
	REAL u(m,n),bb(m),resold(:),resnew(:),
     *     x(NPMAX),xold(:),xnew(:),max,sumresold,sumresnew,
     *     sumbb

	ALLOCATABLE resold,resnew,xold,xnew
	ALLOCATE (resold(m),resnew(m),xold(NPMAX),xnew(NPMAX))

	 write(*,'('' m,n='',2i4)')m,n

	xold(1)=0.
	do 1 j=1,m
	 xold(1)=xold(1)+bb(j)/n/m
1     continue
	if(n.EQ.1)then
	 xnew(1)=xold(1)
	 goto 8
	endif

	do 2 i=2,n
	 xold(i)=xold(1)
2     continue
	do 3 j=1,m
	 resold(j)=0.
	 do 4 i=1,n
	   resold(j)=resold(j)+u(j,i)*xold(i)
4      continue
	 resold(j)=resold(j)-bb(j)
3     continue

	goto 10

	 write(*,'('' x'')')

9     do 11 i=1,n
	 xold(i)=xnew(i)
	 write(*,*)xold(i),xnew(i)
11    continue
	 write(*,'('' res'')')
	do 12 j=1,m
	 resold(j)=resnew(j)
	 write(*,*)resold(j),resnew(j)
12    continue
	 pause

10    max=0.
	do 5 j=1,m
	 if(ABS(resold(j)).GE.max)then
	   max=ABS(resold(j))
	   p=j
	 endif
5     continue

	pp=p

c      if(p.GT.n)pp=MOD(p,n)+1
	if(p.GT.n)then
	 phelp=RND()*n
	write(*,*)n,phelp
	 pp=phelp+1
	endif
	 write(*,'('' p,pp='',2i4)')p,pp

	do 6 i=1,n
	 if(i.EQ.pp)then
	   xnew(i)=xold(i)-resold(p)/u(p,pp)
	   if(xnew(i).LT.0)xnew(i)=0.
	 else
	   xnew(i)=xold(i)
	 endif
6     continue
	sumresold=0.
	sumresnew=0.
	sumbb=0.
	do 7 j=1,m
	   resnew(j)=resold(j)-u(j,pp)*resold(p)/u(p,pp)
	   sumresold=sumresold+ABS(resold(j))
	   sumresnew=sumresnew+ABS(resnew(j))
	   sumbb=sumbb+ABS(bb(j))
7     continue
	if(sumresnew/sumresold.LT.1.005.AND.sumresnew/sumbb.LT.0.1)then
	  goto 8
	else
	  goto 9
	endif

8     do 13 i=1,n
	 x(i)=xnew(i)
13    continue
	DEALLOCATE (resold,resnew,xold,xnew)
	RETURN
	END
C*****************************************************************
	SUBROUTINE svbksb(u,w,v,m,n,mp,np,b,x,
     *NDMAX,MAMAX,NCMAX,MPMAX,NPMAX)

      INTEGER m,mp,n,np,NMAX
	REAL b(mp),u(mp,np),v(np,np),w(np),x(NPMAX)
	PARAMETER (NMAX=1000)
      INTEGER i,j,jj
      REAL s,tmp(NMAX)
      do 12 j=1,n
        s=0.
        if(w(j).ne.0.)then
          do 11 i=1,m
            s=s+u(i,j)*b(i)
11        continue
          s=s/w(j)
        endif
        tmp(j)=s
12    continue
      do 14 j=1,n
        s=0.
        do 13 jj=1,n
          s=s+v(j,jj)*tmp(jj)
13      continue
        x(j)=s
14    continue
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software Y"1i41)jD#ks<,1'i0S%.
C**********************************************************************
	SUBROUTINE svdcmp(a,m,n,mp,np,w,v)
      INTEGER m,mp,n,np,NMAX
      REAL a(mp,np),v(np,np),w(np)
	PARAMETER (NMAX=1000)
CU    USES pythag
      INTEGER i,its,j,jj,k,l,nm
      REAL anorm,c,f,g,h,s,scale,x,y,z,rv1(NMAX),pythag
      g=0.0
      scale=0.0
      anorm=0.0
      do 25 i=1,n
        l=i+1
        rv1(i)=scale*g
        g=0.0
        s=0.0
        scale=0.0
        if(i.le.m)then
          do 11 k=i,m
            scale=scale+abs(a(k,i))
11        continue
          if(scale.ne.0.0)then
            do 12 k=i,m
              a(k,i)=a(k,i)/scale
              s=s+a(k,i)*a(k,i)
12          continue
            f=a(i,i)
            g=-sign(sqrt(s),f)
            h=f*g-s
            a(i,i)=f-g
            do 15 j=l,n
              s=0.0
              do 13 k=i,m
                s=s+a(k,i)*a(k,j)
13            continue
              f=s/h
              do 14 k=i,m
                a(k,j)=a(k,j)+f*a(k,i)
14            continue
15          continue
            do 16 k=i,m
              a(k,i)=scale*a(k,i)
16          continue
          endif
        endif
        w(i)=scale *g
        g=0.0
        s=0.0
        scale=0.0
        if((i.le.m).and.(i.ne.n))then
          do 17 k=l,n
            scale=scale+abs(a(i,k))
17        continue
          if(scale.ne.0.0)then
            do 18 k=l,n
              a(i,k)=a(i,k)/scale
              s=s+a(i,k)*a(i,k)
18          continue
            f=a(i,l)
            g=-sign(sqrt(s),f)
            h=f*g-s
            a(i,l)=f-g
            do 19 k=l,n
              rv1(k)=a(i,k)/h
19          continue
            do 23 j=l,m
              s=0.0
              do 21 k=l,n
                s=s+a(j,k)*a(i,k)
21            continue
              do 22 k=l,n
                a(j,k)=a(j,k)+s*rv1(k)
22            continue
23          continue
            do 24 k=l,n
              a(i,k)=scale*a(i,k)
24          continue
          endif
        endif
        anorm=max(anorm,(abs(w(i))+abs(rv1(i))))
25    continue
      do 32 i=n,1,-1
        if(i.lt.n)then
          if(g.ne.0.0)then
            do 26 j=l,n
              v(j,i)=(a(i,j)/a(i,l))/g
26          continue
            do 29 j=l,n
              s=0.0
              do 27 k=l,n
                s=s+a(i,k)*v(k,j)
27            continue
              do 28 k=l,n
                v(k,j)=v(k,j)+s*v(k,i)
28            continue
29          continue
          endif
          do 31 j=l,n
            v(i,j)=0.0
            v(j,i)=0.0
31        continue
        endif
        v(i,i)=1.0
        g=rv1(i)
        l=i
32    continue
      do 39 i=min(m,n),1,-1
        l=i+1
        g=w(i)
        do 33 j=l,n
          a(i,j)=0.0
33      continue
        if(g.ne.0.0)then
          g=1.0/g
          do 36 j=l,n
            s=0.0
            do 34 k=l,m
              s=s+a(k,i)*a(k,j)
34          continue
            f=(s/a(i,i))*g
            do 35 k=i,m
              a(k,j)=a(k,j)+f*a(k,i)
35          continue
36        continue
          do 37 j=i,m
            a(j,i)=a(j,i)*g
37        continue
        else
          do 38 j= i,m
            a(j,i)=0.0
38        continue
        endif
        a(i,i)=a(i,i)+1.0
39    continue
      do 49 k=n,1,-1
        do 48 its=1,30
          do 41 l=k,1,-1
            nm=l-1
            if((abs(rv1(l))+anorm).eq.anorm)  goto 2
            if((abs(w(nm))+anorm).eq.anorm)  goto 1
41        continue
1         c=0.0
          s=1.0
          do 43 i=l,k
            f=s*rv1(i)
            rv1(i)=c*rv1(i)
            if((abs(f)+anorm).eq.anorm) goto 2
            g=w(i)
            h=pythag(f,g)
            w(i)=h
            h=1.0/h
            c= (g*h)
            s=-(f*h)
            do 42 j=1,m
              y=a(j,nm)
              z=a(j,i)
              a(j,nm)=(y*c)+(z*s)
              a(j,i)=-(y*s)+(z*c)
42          continue
43        continue
2         z=w(k)
          if(l.eq.k)then
            if(z.lt.0.0)then
              w(k)=-z
              do 44 j=1,n
                v(j,k)=-v(j,k)
44            continue
            endif
            goto 3
          endif
          if(its.eq.30) pause 'no convergence in svdcmp'
          x=w(l)
          nm=k-1
          y=w(nm)
          g=rv1(nm)
          h=rv1(k)
          f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y)
          g=pythag(f,1.0)
          f=((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x
          c=1.0
          s=1.0
          do 47 j=l,nm
            i=j+1
            g=rv1(i)
            y=w(i)
            h=s*g
            g=c*g
            z=pythag(f,h)
            rv1(j)=z
            c=f/z
            s=h/z
            f= (x*c)+(g*s)
            g=-(x*s)+(g*c)
            h=y*s
            y=y*c
            do 45 jj=1,n
              x=v(jj,j)
              z=v(jj,i)
              v(jj,j)= (x*c)+(z*s)
              v(jj,i)=-(x*s)+(z*c)
45          continue
            z=pythag(f,h)
            w(j)=z
            if(z.ne.0.0)then
              z=1.0/z
              c=f*z
              s=h*z
            endif
            f= (c*g)+(s*y)
            x=-(s*g)+(c*y)
            do 46 jj=1,m
              y=a(jj,j)
              z=a(jj,i)
              a(jj,j)= (y*c)+(z*s)
              a(jj,i)=-(y*s)+(z*c)
46          continue
47        continue
          rv1(l)=0.0
          rv1(k)=f
          w(k)=x
48      continue
3       continue
49    continue
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software Y"1i41)jD#ks<,1'i0S%.
c*************************************************************************
	FUNCTION pythag(a,b)
      REAL a,b,pythag
      REAL absa,absb
	absa=abs(a)
      absb=abs(b)
      if(absa.gt.absb)then
        pythag=absa*sqrt(1.+(absb/absa)**2)
      else
        if(absb.eq.0.)then
          pythag=0.
        else
          pythag=absb*sqrt(1.+(absa/absb)**2)
        endif
      endif
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software Y"1i41)jD#ks<,1'i0S%.
C --------------------------------------------------------------------------
C ----NEW TEXTURE CORRECTION FOR INCLINED TEXTURE AND ASYMMETRIC DIFFRACTION
C --------------------------------------------------------------------------
	FUNCTION TEXTUR(ph,a,b,rn,RG,is)
C	***************************************************************
C	Function calculates the new texture correction formula for
C     inclined textures.
C	Dr. Radovan Cerny, Charles University,     Prague, 9.7.1992
C	***************************************************************
	REAL THETA,bt,G,n,phiHh,alpha0,P(721),PD(721)
	INTEGER switch,points,step,swi
	COMMON /data/ THETA,bt,G,n,phiHh,alpha0,points,step
	COMMON /PHKLdata/ P,PD
	COMMON /sw/ switch,swi
	EXTERNAL SCANPHKL,RNORM
	dvepi=6.2831853
	IF(bt.NE.b.OR.G.NE.RG.OR.n.NE.rn.OR.(swi.NE.2.AND.is.EQ.2))
     *       THEN
C		write(*,100) is,swi,switch,G
C100		format(1x,'is=',i3,' swi=',i3,' switch=',i3,' G=',f6.2)
		bt=b
		n=rn
		G=RG
		switch=is		
		IF(bt.LE.2.OR.G.GT.20.OR.n.GT.5) THEN
			points=720
			step=4
		ELSE
			points=180
			step=1
		ENDIF
		CALL PHKLROT
	ENDIF
	phiHh=ph
	alpha0=a
	switch=is		
C	write(*,101) phiHh,is,swi,switch,G
C101	format(' 2.pass ','phiHh= ',f5.2,'is=',i3,' swi=',i3,
C     .  ' switch=',i3,' G=',f6.2)
	IF (SWITCH.EQ.1) THEN
	    TEXTUR=GAUSSINT(SCANPHKL,0.,dvepi)/dvepi
     *          /GAUSSINT(RNORM,0.,dvepi/4)
	ELSE
	    SCANP1=GAUSSINT(SCANPHKL,0.,dvepi)/dvepi
	    RNOR2=GAUSSINT(RNORM,0.,dvepi/4)
	    SWITCH=1
	    SCANP2=GAUSSINT(SCANPHKL,0.,dvepi)/dvepi
	    RNOR1=GAUSSINT(RNORM,0.,dvepi/4)
	    SWITCH=2
	    TEXTUR=(RNOR1*SCANP1+RNOR2*SCANP2)/RNOR1**2
	END IF
c      TEXTUR=0.12+0.88*TEXTUR
	RETURN
	END
C     **************************************************************
	FUNCTION RNORM(x)
	REAL THETA,bt,G,n,phiHh,alpha0,P(721),PD(721)
	INTEGER switch,points,step,swi
	COMMON /data/ THETA,bt,G,n,phiHh,alpha0,points,step
	COMMON /sw/ switch,swi
	IF(switch.EQ.1)THEN
	RNORM = EXP(-G*(1-(COS(x))**n))*SIN(x)
	ELSE
	RNORM = (1-(COS(x))**n)*EXP(-G*(1-(COS(x))**n))*SIN(x)
	END IF
	RETURN
	END
C     ***************************************************************
	FUNCTION GAUSSINT(F,a,b)
C     ***************************************************************
C     Gauss 180-point integration of function F(x).
C	***************************************************************
	REAL F,a,b,y1,y2,gauss,x(90),w(90)
	INTEGER i
	DATA (  x(i)                    ,           w(i)     , i=1,90)/
C     1      0.048307665687738316235       ,     0.096540088514727800567,
C     2      0.144471961582796493485       ,     0.095638720079274859419,
C     3      0.239287362252137074545       ,     0.093844399080804565639,
C     4      0.331868602282127649780       ,     0.091173878695763884713,
C     5      0.421351276130635345364       ,     0.087652093004403811143,
C     6      0.506899908932229390024       ,     0.083311924226946755222,
C     7      0.587715757240762329041       ,     0.078193895787070306472,
C     8      0.663044266930215200975       ,     0.072345794108848506225,
C     9      0.732182118740289680387       ,     0.065822222776361846838,
C     A      0.794483795967942406963       ,     0.058684093478535547145,
C     B      0.849367613732569970134       ,     0.050998059262376176196,
C     C      0.896321155766052123965       ,     0.042835898022226680657,
C     D      0.934906075937739689171       ,     0.034273862913021433103,
C     E      0.964762255587506430774       ,     0.025392065309262059456,
C     F      0.985611511545268335400       ,     0.016274394730905670605,
C     G      0.997263861849481563545       ,     0.007018610009470096600/
     C         8.702106E-03    ,    1.740412E-02    ,
     C         2.610422E-02    ,    1.739907E-02    ,
     C         4.349834E-02    ,    1.738860E-02    ,
     C         6.087928E-02    ,    1.737239E-02    ,
     C         7.824179E-02    ,    1.735146E-02    ,
     C         9.558066E-02    ,    1.732528E-02    ,
     C         1.128904E-01    ,    1.729336E-02    ,
     C         1.301661E-01    ,    1.725658E-02    ,
     C         1.474023E-01    ,    1.721444E-02    ,
     C         1.645939E-01    ,    1.716805E-02    ,
     C         1.817355E-01    ,    1.711468E-02    ,
     C         1.988222E-01    ,    1.705734E-02    ,
     C         2.158488E-01    ,    1.699447E-02    ,
     C         2.328096E-01    ,    1.692654E-02    ,
     C         2.497002E-01    ,    1.685375E-02    ,
     C         2.665153E-01    ,    1.677515E-02    ,
     C         2.832493E-01    ,    1.669235E-02    ,
     C         2.998976E-01    ,    1.660390E-02    ,
     C         3.164550E-01    ,    1.651046E-02    ,
     C         3.329166E-01    ,    1.641238E-02    ,
     C         3.492773E-01    ,    1.630837E-02    ,
     C         3.655323E-01    ,    1.620048E-02    ,
     C         3.816764E-01    ,    1.608745E-02    ,
     C         3.977051E-01    ,    1.596901E-02    ,
     C         4.136132E-01    ,    1.584621E-02    ,
     C         4.293963E-01    ,    1.571843E-02    ,
     C         4.450491E-01    ,    1.558621E-02    ,
     C         4.605671E-01    ,    1.544901E-02    ,
     C         4.759455E-01    ,    1.530740E-02    ,
     C         4.911797E-01    ,    1.516044E-02    ,
     C         5.062651E-01    ,    1.501033E-02    ,
     C         5.211973E-01    ,    1.485342E-02    ,
     C         5.359716E-01    ,    1.469385E-02    ,
     C         5.505835E-01    ,    1.452900E-02    ,
     C         5.650285E-01    ,    1.436074E-02    ,
     C         5.793026E-01    ,    1.418679E-02    ,
     C         5.934010E-01    ,    1.400965E-02    ,
     C         6.073197E-01    ,    1.382748E-02    ,
     C         6.210544E-01    ,    1.364182E-02    ,
     C         6.346009E-01    ,    1.345079E-02    ,
     C         6.479554E-01    ,    1.325609E-02    ,
     C         6.611134E-01    ,    1.305870E-02    ,
     C         6.740714E-01    ,    1.285634E-02    ,
     C         6.868250E-01    ,    1.264989E-02    ,
     C         6.993706E-01    ,    1.244071E-02    ,
     C         7.117042E-01    ,    1.222674E-02    ,
     C         7.238224E-01    ,    1.200916E-02    ,
     C         7.357212E-01    ,    1.178754E-02    ,
     C         7.473972E-01    ,    1.156372E-02    ,
     C         7.588466E-01    ,    1.133473E-02    ,
     C         7.700664E-01    ,    1.110408E-02    ,
     C         7.810529E-01    ,    1.086846E-02    ,
     C         7.918025E-01    ,    1.063107E-02    ,
     C         8.023127E-01    ,    1.038883E-02    ,
     C         8.125796E-01    ,    1.014397E-02    ,
     C         8.226005E-01    ,    9.896754E-03    ,
     C         8.323719E-01    ,    9.646029E-03    ,
     C         8.418913E-01    ,    9.392838E-03    ,
     C         8.511558E-01    ,    9.135622E-03    ,
     C         8.601624E-01    ,    8.877348E-03    ,
     C         8.689083E-01    ,    8.615205E-03    ,
     C         8.773911E-01    ,    8.350058E-03    ,
     C         8.856084E-01    ,    8.082921E-03    ,
     C         8.935571E-01    ,    7.813577E-03    ,
     C         9.012350E-01    ,    7.542177E-03    ,
     C         9.086400E-01    ,    7.267451E-03    ,
     C         9.157699E-01    ,    6.991300E-03    ,
     C         9.226225E-01    ,    6.713357E-03    ,
     C         9.291950E-01    ,    6.432091E-03    ,
     C         9.354867E-01    ,    6.150329E-03    ,
     C         9.414949E-01    ,    5.865652E-03    ,
     C         9.472176E-01    ,    5.579985E-03    ,
     C         9.526535E-01    ,    5.292066E-03    ,
     C         9.578009E-01    ,    5.002728E-03    ,
     C         9.626579E-01    ,    4.711291E-03    ,
     C         9.672237E-01    ,    4.419836E-03    ,
     C         9.714963E-01    ,    4.125773E-03    ,
     C         9.754745E-01    ,    3.830137E-03    ,
     C         9.791573E-01    ,    3.534990E-03    ,
     C         9.825435E-01    ,    3.237800E-03    ,
     C         9.856320E-01    ,    2.939536E-03    ,
     C         9.884220E-01    ,    2.639816E-03    ,
     C         9.909126E-01    ,    2.340911E-03    ,
     C         9.931031E-01    ,    2.039973E-03    ,
     C         9.949927E-01    ,    1.738935E-03    ,
     C         9.965807E-01    ,    1.436925E-03    ,
     C         9.978669E-01    ,    1.135092E-03    ,
     C         9.988511E-01    ,    8.323371E-04    ,
     C         9.995325E-01    ,    5.306011E-04    ,
     C         9.999114E-01    ,    2.278413E-04/
C       2.435030E-02  ,    4.869149E-02  ,
C       7.299300E-02  ,    4.857572E-02  ,
C       1.214627E-01  ,    4.834482E-02  ,
C       1.696445E-01  ,    4.799954E-02  ,
C       2.174237E-01  ,    4.754064E-02  ,
C       2.646871E-01  ,    4.696765E-02  ,
C       3.113229E-01  ,    4.628531E-02  ,
C       3.572200E-01  ,    4.549139E-02  ,
C       4.022702E-01  ,    4.459032E-02  ,
C       4.463658E-01  ,    4.358399E-02  ,
C       4.894033E-01  ,    4.247323E-02  ,
C       5.312794E-01  ,    4.126259E-02  ,
C       5.718957E-01  ,    3.995387E-02  ,
C       6.111553E-01  ,    3.855056E-02  ,
C       6.489656E-01  ,    3.705504E-02  ,
C       6.852361E-01  ,    3.547222E-02  ,
C       7.198820E-01  ,    3.380480E-02  ,
C       7.528199E-01  ,    3.205802E-02  ,
C       7.839724E-01  ,    3.023459E-02  ,
C       8.132654E-01  ,    2.833976E-02  ,
C       8.406293E-01  ,    2.637777E-02  ,
C       8.659993E-01  ,    2.435239E-02  ,
C       8.893153E-01  ,    2.227057E-02  ,
C       9.105220E-01  ,    2.013475E-02  ,
C       9.295692E-01  ,    1.795184E-02  ,
C       9.464113E-01  ,    1.572618E-02  ,
C       9.610087E-01  ,    1.346273E-02  ,
C       9.733267E-01  ,    1.116824E-02  ,
C       9.833360E-01  ,    8.846750E-03  ,
C       9.910134E-01  ,    6.504408E-03  ,
C       9.963400E-01  ,    4.146992E-03  ,
C       9.993050E-01  ,    1.783196E-03/
	GAUSS=0
	DO 1 i=1,90
	y1=(b-a)/2*x(i)+(b+a)/2
	y2=-(b-a)/2*x(i)+(b+a)/2
	GAUSS=GAUSS+(F(y1)+F(y2))*w(i)
1     CONTINUE
	GAUSSINT=(b-a)/2*GAUSS
	RETURN
	END
C     ****************************************************************
	FUNCTION PHKL(x)
	REAL THETA,bt,G,n,help,x
	INTEGER switch,points,step,swi
	COMMON /data/ THETA,bt,G,n,phiHh,alpha0,points,step
	COMMON /sw/ switch,swi
	dvepi=6.2831853
C     WRITE(*,100) bt,THETA,x
C100  FORMAT(1X,'BETA= ',F5.2,' THETA= ',F5.2,' x= ',f5.2)
	a=COS(bt/360*dvepi)*COS(THETA/360*dvepi)
	b=SIN(bt/360*dvepi)*SIN(THETA/360*dvepi)
	IF (switch.EQ.2) THEN
		help=(ABS(a-b*SIN(x)))**n-1
	ELSE
		help=1.
	ENDIF
	PHKL=EXP(-G)/dvepi*EXP(G*(ABS(a-b*SIN(x)))**n)*help
	RETURN
	END
C     *****************************************************************
	SUBROUTINE PHKLROT
	REAL THETA,bt,G,n,phiHh,alpha0,P(721),PD(721)
	INTEGER switch,points,step,swi
	COMMON /data/ THETA,bt,G,n,phiHh,alpha0,points,step
	COMMON /PHKLdata/ P,PD
	COMMON /sw/ switch,swi
	EXTERNAL PHKL
	dvepi=6.2831853
	IF(bt.EQ.0.) GOTO 2
	DO 1 i=0,points
	THETA=FLOAT(i)/step
	IF(switch.EQ.2) THEN
		swi=2
		PD(i+1)=GAUSSINT(PHKL,0.,dvepi)
	ELSE
		swi=1
		P(i+1)=GAUSSINT(PHKL,0.,dvepi)
	ENDIF
1	CONTINUE
2	RETURN
	END
C	*****************************************************************
	FUNCTION SCANPHKL(x)
	REAL THETA,bt,G,n,phiHh,alpha0,x,P(721),PD(721)
	INTEGER switch,points,step,swi
	COMMON /data/ THETA,bt,G,n,phiHh,alpha0,points,step
	COMMON /PHKLdata/ P,PD
	COMMON /sw/ switch,swi
	dvepi=6.2831853
C	write(*,100) x,phiHh
C100	format(1x,' x= ',f5.2,' phiHh= ',f5.2)
	THETA=-SIN(x)*SIN(phiHh)*SIN(alpha0)+
     1	      COS(phiHh)*COS(alpha0)
	THETA=ACOS(THETA)/dvepi*360
	IF(bt.EQ.0.) THEN
	   SCANPHKL=dvepi*PHKL(0.)
	ELSE
	   DO 1 i=0,points
	   IF(THETA.LE.FLOAT(i)/step) GOTO 2
1	   CONTINUE
2	   IF(switch.EQ.2) THEN
		    SCANPHKL=(PD(i+1)-PD(i))*step/(i-(i-1))*
     1              (THETA-FLOAT(i-1)/step)+PD(i)
	   	    IF(THETA.EQ.0.) SCANPHKL=PD(1)
	   ELSE	  
		    SCANPHKL=(P(i+1)-P(i))*step/(i-(i-1))*
     1              (THETA-FLOAT(i-1)/step)+P(i)
	   	    IF(THETA.EQ.0.) SCANPHKL=P(1)
	   ENDIF
	ENDIF
	RETURN
	END

*     ************************
	SUBROUTINE texsphh(dc,hh,kk,ll,
     *          jump,geom,texturecor,coef,texder,norder)
*     ************************
*     Texture correction with symmetrized spherical harmonics.
*
*     This subroutine is a modification of the code from Helsinki.
*     Modified by
*                   Radovan Cerny
*                   Laboratoire de Cristalographie
*                   University of Geneva
*     6.11.2000
*     **************************************************************
*     Input: MILERIND(3) - hkl
*            DC(6)       - a,b,c,alpha,beta,gamma; angles in deg
*            JUMP        - Laue symmetry:
*                           JUMP                    Laue group
*                            1                         -1
*                            2                         2/m
*                            3                         mmm a-axis
*                            4                         mmm b-axis
*                            5                         mmm c-axis
*                            6                         4/m
*                            7                         4/mmm
*                            8                         -3   in hex. lattice
*                            9                         -3m1 in hex. lattice
*                                                      -31m in hex. lattice
*                                                      -3m  in hex. lattice
*                            10                        6/m
*                            11                        6/mmm
*                            12                        m-3
*                            13                        m-3m
*            GEOM        - Cosinus od the angle between
*                          the texture axis (spinning axis)
*                          and the diffraction vector :
*                           GEOM                       Geometry
*                             1       Symmetric reflection, no chi-inclination
*                           COS(chi)  Symmetric reflection, chi-inclination
*                             0       Transmission or Debye-Scherrer in
*                                     eqatorial plane
*                   COS(theta-gamma)  Asymetric, no chi-inclination (Guinier):
*                                       reflection : gamma = angle between
*                                     a plane perpendicular to the texture
*                                     (spinning) axis (usually specimen surface)
*                                     and incident beam
*                                       transmission : gamma = - gamma
*           COS(theta-gamma)*COS(chi) Asymetric, chi-inclination
*                                     Image plate data can be described by this
*                                     option as transmission data with :
*                                       gamma = - gamma
*                                       chi = azimuthal angle on image plate,
*                                             zero corresponds to the direction
*                                             of the texture (spinning) axis
*                                             projection on the image plate.
*     ***********************************************************************

	INTEGER hh,kk,ll,jump,norder
	REAL dc(6),geom,texturecor,coef(15),texder(15)

	CALL coord(hh,kk,ll,dc,jump,costhe,cosphi)

	CALL textur2(costhe,cosphi,geom,jump,texturecor,coef,texder,
     *             norder)

	RETURN
	END
*     ***********************
	SUBROUTINE COORD(hh,kk,ll,dc,jump,costhe,cosphi)
*     ***********************
*
*     This subroutine sets the correct cartesian coordinate orientation
*     according to space symmetry.
*
	INTEGER hh,kk,ll
	INTEGER JUMP
	REAL H,K,L,pi,DC(6)
      REAL A,B,C,COSALF,COSBET,COSGAM,SINALF,SINBET,SINGAM
      REAL COSPHI,COSTHE

	pi=3.141592654

	H=FLOAT(hh)
	K=FLOAT(kk)
	L=FLOAT(ll)
      A=DC(1)
      B=DC(2)
      C=DC(3)
	COSALF=COS(DC(4)/180.*pi)
	COSBET=COS(DC(5)/180.*pi)
	COSGAM=COS(DC(6)/180.*pi)
      SINALF=SQRT(1-COSALF**2)
      SINBET=SQRT(1-COSBET**2)
      SINGAM=SQRT(1-COSGAM**2)
*
*     Branch according to the structure
*
      GO TO (10,20,30,30,30,60,60,80,80,80,80,120,120),JUMP
      GO TO 150
*
*     Triclinic         -- Space groups 1   ...   2
*
  10  S11=(B*C*SINALF)**2
      S12=A*B*C*C*(COSALF*COSBET-COSGAM)
      S22=(A*C*SINBET)**2
      S33=(A*B*SINGAM)**2
      S23=A*A*B*C*(COSBET*COSGAM-COSALF)
      S13=A*B*B*C*(COSGAM*COSALF-COSBET)
*
      D=SQRT(S33*(S11*H*H+S22*K*K+S33*L*L+2*(S23*K*L+S13*H*L+S12*H*K)))
      CT=(S33*L+S23*K+S13*H)/D
      D=SQRT(S11*(S11*H*H+S22*K*K+2*S12*H*K))
      CP=(S11*H+S12*K)/D
      GO TO 140
*
*     Monoclinic        -- Space groups 3   ...  15  -  b-axis unique
*
  20  D=(H/A)**2+(L/C)**2-2*H*L*COSBET/(A*C)
      IF(D.EQ.0) THEN
         CP=1
      ELSE
         CP=(H/A-L*COSBET/C)/SQRT(D)
      END IF
      CT=K*SINBET/(B*SQRT(D+(K*SINBET/B)**2))
      GO TO 140
*
*     Orthorhombic      -- Space groups 16  ...  74
*
  30  D=SQRT((H/A)**2+(K/B)**2+(L/C)**2)
      IF(JUMP.EQ.3) THEN
*                                  Case A: a-axis unique
         CT=H/(A*D)
         OVER=K/B
         UNDER=SQRT((K/B)**2+(L/C)**2)
      ELSE IF(JUMP.EQ.4) THEN
*                                  Case B: b-axis unique
         CT=K/(B*D)
         OVER=L/C
         UNDER=SQRT((H/A)**2+(L/C)**2)
*                                  Case C: c-axis unique
      ELSE
         CT=L/(C*D)
         OVER=H/A
         UNDER=SQRT((H/A)**2+(K/B)**2)
      END IF
      IF(UNDER.EQ.0) THEN
         CP=1
      ELSE
         CP=OVER/UNDER
      END IF
      GO TO 140
*
*     Tetragonal          -- Space groups 75  ... 142
*
   60 D=SQRT((H*H+K*K)/(A*A)+(L/C)**2)
      CT=L/(C*D)
      IF(H*H+K*K.EQ.0) THEN
         CP=1
      ELSE
         CP=H/SQRT(H*H+K*K)
      END IF
      GO TO 140
*
*     Trigonal            -- Space groups 143 ... 167
*     Hexagonal           -- Space groups 168 ... 194
*
   80 D=H*H+K*K+H*K
      CT=0.8660254038*A*L/(C*SQRT(D+3*(A*L/(2*C))**2))
      IF(D.EQ.0) THEN
         CP=1
      ELSE
         CP=(H+0.5*K)/SQRT(D)
      END IF
      GO TO 140
*
*     Cubic               -- Space groups 195 ... 230
*
  120 D=H*H+K*K
      CT=L/SQRT(D+L*L)
      IF(D.EQ.0) THEN
         CP=1
      ELSE
         CP=H/SQRT(D)
      END IF
*
*     Get rid of rounding errors.
*
  140 IF(ABS(CT).GT.1) CT=SIGN(1.0,CT)
      IF(ABS(CP).GT.1) CP=SIGN(1.0,CP)
	COSTHE=CT
	COSPHI=CP
      RETURN
  150 WRITE(6,1000) JUMP
      STOP
 1000 FORMAT(' In COORD: Invalid symmetry; JUMP =',I5,'.')
	END

*     **************************
	SUBROUTINE TEXTUR2(costhe,cosphi,geom,jump,texturecor,coef,texder,
     *                   norder)
*     **************************
*
	INTEGER ORDER(15,13),norder
	REAL COSPHI,COSTHE,COEF(15),TEXDER(15),texturecor

*   JUMP               Indices of Y
      DATA ((order(i,j),i=1,15),j=1,3)/
     1      20, 21,-21, 22,-22, 40, 41,-41, 42,-42, 43,-43, 44,-44, 60,
     2      20, 22,-22, 40, 42,-42, 44,-44, 60, 62,-62, 64,-64, 66,-66,
     3      20, 22, 40, 42, 44, 60, 62, 64, 66, 80,100,120,  0,  0,  0/
      DATA ((order(i,j),i=1,15),j=4,6)/
     4      20, 22, 40, 42, 44, 60, 62, 64, 66, 80,100,120,  0,  0,  0,
     5      20, 22, 40, 42, 44, 60, 62, 64, 66, 80,100,120,  0,  0,  0,
     6      20, 40, 44,-44, 60, 64,-64, 80,100,120,  0,  0,  0,  0,  0/
      DATA ((order(i,j),i=1,15),j=7,9)/
     7      20, 40, 44, 60, 64, 80,100,120,  0,  0,  0,  0,  0,  0,  0,
     8      20, 40, 43,-43, 60, 63,-63, 66,-66, 80,100,120,  0,  0,  0,
     9      20, 40, 43, 60, 63, 66, 80,100,120,  0,  0,  0,  0,  0,  0/
      DATA ((order(i,j),i=1,15),j=10,12)/
     A      20, 40, 60, 66,-66, 80,100,120,  0,  0,  0,  0,  0,  0,  0,
     B      20, 40, 60, 66, 80,100,120,  0,  0,  0,  0,  0,  0,  0,  0,
     C      40, 61, 62, 81,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0/
      DATA (order(i,13),i=1,15)/
     D      40, 61, 81,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0/
*
      DO 10 I=1,15
*
*     Pick up indices lms for spherical harmonic functions
*     and branch out when l=0 is encountered.
*
      L=ABS(ORDER(I,JUMP)/10)
	IF(L .EQ. 0)GO TO 20
      M=MOD(ORDER(I,JUMP),10)
*                             Spherical harmonic functions come from Y
      IF(JUMP .LT. 12) THEN
	TEXDER(I)=Y(L,0,GEOM,0.0)*Y(L,M,COSTHE,COSPHI)
*
*                             Cubic harmonic functions come from CUBHAR
      ELSE
	TEXDER(I)=Y(L,0,GEOM,0.0)*CUBHAR(L,M,COSTHE,COSPHI)
      END IF
   10 CONTINUE
*
*     Texture correction coefficient is sum of products
*     parameter * harmonic function
*
   20 TEX=1
	DO 30 M=1,I-1
      TEX=TEX+TEXDER(M)*COEF(M)
   30 CONTINUE
	texturecor=TEX
      RETURN
      END
*     *****************************
      FUNCTION Y(L,N,COSTHE,COSPHI)
*     *****************************
*
      INTEGER L,N
      REAL COSTHE,COSPHI,Y
*
*     This subroutine evaluates spherical harmonic functions
*
*     Y   (theta,phi),   where m=ABS(N) and s is sign of N.
*      lms
*
      INTEGER M
      REAL SINTHE
*
      M=ABS(N)
      SINTHE=SQRT(1-COSTHE**2)
*
*     First operate with l only - branch according to it.
*     Odd values of l are invalid for powder neutron diffraction.
*
      GO TO (900,20,900,40,900,60,900,80,900,100,900,120,900),L
      GO TO 900
*                             l = 2
   20 IF(M .EQ. 0) THEN
         YL=(3*COSTHE**2-1)/2
      ELSE IF(M .EQ. 1) THEN
         YL=2*COSTHE*SINTHE
      ELSE
         YL=SINTHE**2
      END IF
      GO TO 200
*                             l = 4
   40 IF(M .EQ. 0) THEN
         YL=(COSTHE**2*(35*COSTHE**2-30)+3)/8
      ELSE IF(M. EQ. 1) THEN
         YL=0.94695*(7*COSTHE**2-3)*COSTHE*SINTHE
      ELSE IF(M .EQ. 2) THEN
         YL=7*(7*COSTHE**2-1)*SINTHE**2/9
      ELSE IF(M .EQ. 3) THEN
         YL=3.079201436*COSTHE*SINTHE**3
      ELSE
         YL=SINTHE**4
      END IF
      GO TO 200
*                             l = 6
   60 IF(M .EQ. 0) THEN
         YL=(((231*COSTHE**2-315)*COSTHE**2+105)*COSTHE**2-5)/16
      ELSE IF(M .EQ. 1) THEN
         GO TO 900
      ELSE IF(M .EQ. 2) THEN
         YL=0.64549*((33*COSTHE**2-18)*COSTHE**2+1)*SINTHE**2
      ELSE IF(M .EQ. 3) THEN
         YL=1.41685*(11*COSTHE**2-3)*COSTHE*SINTHE**3
      ELSE IF(M .EQ. 4) THEN
         YL=0.81675*(11*COSTHE**2-1)*SINTHE**4
      ELSE IF(M .EQ. 5) THEN
         GO TO 900
      ELSE
         YL=SINTHE**6
      END IF
      GO TO 200
*                             l = 8
   80 IF(M .EQ. 0) THEN
         YL=((((6435*COSTHE**2-12012)*COSTHE**2+6930)*COSTHE**2-1260)*
     *          COSTHE**2+35)/128
      ELSE
         GO TO 900
      END IF
      GO TO 200
*                             l = 10
  100 IF(M.EQ.0) THEN
	 YL=(((((46189*COSTHE**2-109395)*COSTHE**2+90090)*COSTHE**2
     *      -30030)*COSTHE**2+3465)*COSTHE**2-63)/256
      ELSE
	 GO TO 900
      END IF
      GO TO 200
*                             L = 12
  120 IF(M.EQ.0) THEN
	 YL=((((((676039*COSTHE**2-1939938)*COSTHE**2+2078505)*COSTHE**2
     *      -1021020)*COSTHE**2+225225)*COSTHE**2-18018)*COSTHE**2+231)/
     *      1024
      ELSE
 	 GO TO 900
      END IF
*
*     Now m, that is, phi-dependence is added.
*
  200 IF(N .GT. 0) THEN
         Y=YL*COS(M*ACOS(COSPHI))
      ELSE IF(N .LT. 0) THEN
         Y=YL*SIN(M*ACOS(COSPHI))
      ELSE
         Y=YL
      END IF
*
      RETURN
*
*     Error: L is illegal.
*
  900 WRITE(6,1000) L,M
      STOP
 1000 FORMAT(' In Y: Invalid indices; L =',I3,' M =',I3,'.')
      END
*     **********************************
      FUNCTION CUBHAR(L,M,COSTHE,COSPHI)
*     **********************************
*
      INTEGER L,M
      REAL K,COSTHE,COSPHI
*
*     This subroutine evaluates cubic harmonic functions
*     K  (theta,phi).
*      lm
*
*     Global statements
*
      X2 = (1-COSTHE**2)*COSPHI**2
      X4 = X2*X2
      X6 = X4*X2
      Y2 = (1-COSTHE**2)*(1-COSPHI**2)
      Y4 = Y2*Y2
      Y6 = Y4*Y2
      Z2 = COSTHE**2
      Z4 = Z2*Z2
      Z6 = Z4*Z2
*
*     Branch with l - values 4, 6 and 8 are allowed.
*
      GO TO (120,120,120,40,120,60,120,80),L
      GO TO 120
   40 K=2.5*(X4+Y4+Z4)-1.5
      GO TO 110
   60 IF(M .EQ. 1) THEN
         K=(693*(X6+Y6+Z6)-945*(X4+Y4+Z4)+270)/32
      ELSE IF(M .EQ. 2) THEN
         K=10.39230485*((Y2-Z2)*X4+(Z2-Z2)*Y4+(X2-Y2)*Z4)
      ELSE
         GO TO 120
      END IF
      GO TO 110
   80 IF(M .EQ. 1) THEN
         K=(195*(X6*X2+Y6*Y2+Z6*Z2)-364*(X6+Y6+Z6)+105*(X4+Y4+Z4)-35)/6
      ELSE
         GO TO 120
      END IF
  110 CUBHAR=K
      RETURN
*
*     Error: L or M is illegal.
*
  120 WRITE(6,1000) L,M
      STOP
 1000 FORMAT('In CUBHAR: Invalid indices; L =',I3,' M =',I3,'.')
      END

C     ****************************************************************
	FUNCTION angPhiHh(h,k,l,ihtx,iktx,iltx,ao,bo,co,alpha,beta,gamma)
C**********************************************************************
C     Function calculate the angle between planes (hkl) and (ihtx,iktx,iltx).
C******************************************************************************
	integer h,k,l,ihtx,iktx,iltx
	real xicos

	pi=3.1415926

	alpharad=alpha/180*pi
      betarad=beta/180*pi
      gammarad=gamma/180*pi
      alphacos=((cos(betarad)*cos(gammarad))-cos(alpharad))/
     1          (sin(betarad)*sin(gammarad))
      betacos=((cos(gammarad)*cos(alpharad))-cos(betarad))/
     1          (sin(gammarad)*sin(alpharad))
      gammacos=((cos(alpharad)*cos(betarad))-cos(gammarad))/
     1          (sin(alpharad)*sin(betarad))

      s=(alpharad+betarad+gammarad)/2
      vol=2*ao*bo*co*sqrt(sin(s)*sin(s-alpharad)*sin(s-betarad)*
     1    sin(s-gammarad))
      arez=bo*co*sin(alpharad)/vol
      brez=co*ao*sin(betarad)/vol
      crez=ao*bo*sin(gammarad)/vol
      qtex1=ihtx**2*arez**2+iktx**2*brez**2+iltx**2*crez**2
      qtex2=2*iktx*iltx*brez*crez*alphacos
      qtex3=2*iltx*ihtx*crez*arez*betacos
      qtex4=2*ihtx*iktx*arez*brez*gammacos
      qtex=qtex1+qtex2+qtex3+qtex4

	    q1=h**2*arez**2+k**2*brez**2+l**2*crez**2
	    q2=2*k*l*brez*crez*alphacos
	    q3=  2*l*h*crez*arez*betacos
	    q4=  2*h*k*arez*brez*gammacos
          q=q1+q2+q3+q4
	  xicos1=h*ihtx*arez**2+k*iktx*brez**2+l*iltx*crez**2
	  xicos2=(k*iltx+l*iktx)*brez*crez*alphacos
	  xicos3=(l*ihtx+h*iltx)*crez*arez*betacos
	  xicos4=(h*iktx+k*ihtx)*arez*brez*gammacos
	  xicos=(xicos1+xicos2+xicos3+xicos4)/sqrt(q*qtex)

		 if (xicos.le.1.00000001) then
		 angPhiHh=acos(xicos)*180/pi
		 IF (angPhiHh.GT.90.) angPhiHh=180.-angPhiHh
		 else
		  goto 30
		 end if

30    return
	end
C*****************************************************************
	FUNCTION rlpg(h,k,l,latpar,wavelength,t)
C******************************************************************
C     Function calculates the Lorentz-polarisation-geometrical correction.
C*************************************************************************
	INTEGER h,k,l
	REAL latpar(6),a,b,c,al,be,ga,wavelength

	rmono=0.8023
	pi=3.1415926

	a=latpar(1)
	b=latpar(2)
	c=latpar(3)
	al=latpar(4)/180*pi
	be=latpar(5)/180*pi
	ga=latpar(6)/180*pi

	d=h*h/a/a*SIN(al)*SIN(al)+k*k/b/b*SIN(be)*SIN(be)+
     *l*l/c/c*SIN(ga)*SIN(ga)+2*h*k/a/b*(COS(al)*COS(be)-COS(ga))+
     *2*k*l/b/c*(COS(be)*COS(ga)-COS(al))+
     *2*l*h/c/a*(COS(ga)*COS(al)-COS(be))
	d=d/(1-COS(al)*COS(al)-
     *COS(be)*COS(be)-COS(ga)*COS(ga)+2*COS(al)*COS(be)*COS(ga))
	d=1/SQRT(d)
	t=ASIN(wavelength/2/d)
	rlpg=(1+rmono*COS(2*t)*COS(2*t))/(1+rmono)/(SIN(2*t)*SIN(t))

	RETURN
	END
C***********************************************************************
