*     *****************
      SUBROUTINE PREPAR
*     *****************
******************************************************
*     This program prepares the profile data.        *
******************************************************

      INTEGER OFFSET
      PARAMETER (OFFSET = 2**16)

      CHARACTER TEXT*80,STAR*1
      LOGICAL MARK
      REAL BAK(100),HW(999),PLOR(999),POS(100),TL(999),T(999),HWOUT(3),
     1     TOUT(3),HKL,SJJ,YOBS,WEIGHT
      INTEGER IOUT(3),ICOUT(3),IHOUT(3,3),ICODE,IH(3,999),IHIGH(10),
     1     INDEX(999),LOW(10),LL(999),NH,M(999),MOUT(3),JCODE(999),IORD
      INCLUDE 'BLANK.FOR'
      COMMON/DATA/CUT,NREFL,NH(3),IALFA,ISTAP,NCHAN,ICODE(999),
     *       HKL(3,999),SJJ(999),YOBS(15000),WEIGHT(15000),IORD(15000)
      COMMON/PARAMS/AINSTR(10),ACRYST(22), ATOM(14,30),
     *              LINSTR(10),LCRYST(22),LATOM(14,30),
     *              PINSTR(10),PCRYST(22),PATOM(14,30)
      COMMON/SPACE/JUMP,A,B,C,COSALF,COSBET,COSGAM,SINALF,SINBET,SINGAM

      DATA SQRLN2/0.83255461116/
      DATA STAR/'*'/
*
*     When angles in 1/10000 cycles, following statement should read:
*     DATA CONST/3183.0988618/  i. e. CONST = 10000/3.141592658...
      DATA CONST/11459.155903/
*                                 Now CONST = 36000/3.141592658...
*     Read and print title.
*
      OPEN(UNIT=3,STATUS='OLD',FORM='FORMATTED'
     *           ,READONLY                          ! VAX/VMS
     *                    )
      OPEN(UNIT=7,STATUS='NEW',FORM='FORMATTED')
      WRITE(7,1000)
      READ(3,'(A)') TEXT
      WRITE(7,'(1X,A)') TEXT
*
      CALL READER(3,TEXT,0,POS,3,IOUT,0,.TRUE.)
      CTHM=POS(1)
      POLAR=POS(2)
      CUT=POS(3)
*
*     Set default value for cut.
*
      IF(CUT.LE.0.OR.CUT.GT.100) CUT=5
      WRITE(7,1001) 2*SQRT(SLABDA)
      WRITE(7,1002)
      WRITE(7,1003) (PINSTR(I),I=3,8)
      WRITE(7,1004) (PINSTR(I),I=9,10)
      WRITE(7,1005) CTHM,POLAR,CUT
*
*     Read in background.
*
      I=1
   10 CALL READER(3,TEXT,0,TOUT,2,IOUT,0,.TRUE.)
      POS(I)=TOUT(1)
      BAK(I)=TOUT(2)
      IF(POS(I).GT.-99) THEN
         I=I+1
         GO TO 10
      END IF
      NB=I-1
      WRITE(7,1006)
      WRITE(7,1007) (POS(I),BAK(I),I=1,NB)
      CALL READER(3,TEXT,0,TOUT,0,IHIGH,2,.TRUE.)
      MREFL=IHIGH(1)
      IHIGH(1)=IHIGH(2)
*
*     When MREFL is negative, the lower and upper limits of those parts
*     of the diagram, which have to be excluded from the calculations,
*     have to be read in.
*
      NCU=0
      IF(MREFL .GT. 0) GO TO 30
      I=1
   20 IF(I.EQ.1) THEN
         LOW(I)=-MREFL
         IHIGH(I)=-IHIGH(I)
      ELSE
         CALL READER(3,TEXT,0,TOUT,0,NH,2,.TRUE.)
         LOW(I)=-NH(1)
         IHIGH(I)=-NH(2)
      END IF
      IF(LOW(I).NE.100) THEN
         I=I+1
         GO TO 20
      END IF
      NCU=I-1
      WRITE(7,1008) (LOW(I),IHIGH(I),I=1,NCU)
      CALL READER(3,TEXT,0,TOUT,0,NH,1,.TRUE.)
      MREFL=NH(1)
   30 IF(NCU.EQ.0) WRITE(7,1009)
*
*     Read in denominators of Miller indices.
*
      CALL READER(3,TEXT,0,TOUT,0,NH,3,.TRUE.)
      WRITE(7,1010)
      WRITE(7,1011) (NH,I=1,3)
*
*     Read reflections and calculate for each the Lorentz factor,
*     the position and the halfwidth.
*
      DO 50 I=1,MREFL
      CALL READER(3,TEXT,0,TOUT,0,JCODE(I),1,.TRUE.)
      CALL READER(3,TEXT,0,TOUT,0,IH(1,I),3,.FALSE.)
      CALL READER(3,TEXT,0,TOUT,0,M(I),1,.FALSE.)
      D=PCRYST(1)*IH(1,I)*IH(1,I)/(NH(1)*NH(1))+
     *  PCRYST(2)*IH(2,I)*IH(2,I)/(NH(2)*NH(2))+
     *  PCRYST(3)*IH(3,I)*IH(3,I)/(NH(3)*NH(3))+
     *  PCRYST(4)*IH(2,I)*IH(3,I)/(NH(2)*NH(3))+
     *  PCRYST(5)*IH(1,I)*IH(3,I)/(NH(1)*NH(3))+
     *  PCRYST(6)*IH(1,I)*IH(2,I)/(NH(1)*NH(2))
      D=SLABDA*D
      PLOR(I)=(1-POLAR+POLAR*CTHM*(1-2*D)**2)/(2*D*SQRT(1-D))
      TG=SQRT(D/(1-D))
      TL(I)=ATAN(TG)*CONST+PINSTR(9)+PINSTR(10)/TG
*
*     Complex part of Voigtian argument (=y) is a function of five
*     halfwidth parameters and Bragg angle:
*
      W2G=SQRT(PINSTR(3)/TG**2+TG*(PINSTR(4)*TG+PINSTR(5))+PINSTR(6))
      Y=SQRLN2*(PINSTR(7)*TG+PINSTR(8)/SQRT(1-D))/W2G
      SJJ(I)=W2G*APPROX(Y)/EXPERF(Y)
   50 INDEX(I)=0
*
*     Sort angles in ascending order.
*
      CALL SORT(TL,LL,MREFL)
      DO 60 I=1,MREFL
      J=LL(I)
      T(I)=TL(J)
   60 HW(I)=SJJ(J)
*
*     Read in profile intensities.
*
      CALL READER(3,TEXT,0,TOUT,0,IOUT,3,.TRUE.)
      IALFA=IOUT(1)
      ISTAP=IOUT(2)
      IOMEG=IOUT(3)
      POS(NB+1)=IOMEG+ISTAP
      BAK(NB+1)=BAK(NB)
      MAXR=MREFL
      PI=0.01*ISTAP*SQRT(CONST)
      JJ=1
      I=1
      MINR=1
      IS=IALFA
*
*     Read in one record of 10 intensities.
*
   70 J=JJ-(JJ-1)/10*10
      JJ9=JJ+9
      IF(J.EQ.1) CALL READER(3,TEXT,0,YOBS(JJ),10,IOUT,0,.TRUE.)
      Y=YOBS(JJ)
*
*     Determine whether intensity is in excluded range.
*
      IORD(JJ)=0
      DO 80 J=1,NCU
      IF(IS.GE.LOW(J).AND.IS.LE.IHIGH(J)) Y=0
   80 YOBS(JJ)=Y
      IF(Y.EQ.0) GO TO 110
      MARK=.FALSE.
*
*     Determine what reflections can contribute to intensity i.e. are
*     within CUT*HALFWIDTH from intensity position.
*
      DO 90 J=MINR,MAXR
      WIDTH=CUT*HW(J)
      ST=T(J)
      MIN=ST-WIDTH+0.5
      MAX=ST+WIDTH+0.5
      IF(IS.GE.MIN.AND.IS.LE.MAX) THEN
         IF(.NOT.MARK) THEN
            I=J
            IORD(JJ)=J
            MARK=.TRUE.
         END IF
         INDEX(J)=-1
      ELSE
         IF(MARK) GO TO 100
      END IF
   90 CONTINUE
      J=MAXR+1
  100 IF(MARK) IORD(JJ)=IORD(JJ) + OFFSET*(J-1)
      MINR=I
  110 IS=IS+ISTAP
      JJ=JJ+1
      IF(IS.LE.IOMEG) GO TO 70
*
*     Determination of total number of reflections within range
*     concerned.
*
      CLOSE(UNIT=3,STATUS='KEEP')
      DO 120 I=1,JJ-1
  120 IF(MOD(IORD(I),OFFSET) .EQ. 0) YOBS(I)=0
      NREFL=0
      DO 130 I=1,MREFL
      IF(INDEX(I).EQ.-1) NREFL=NREFL+1
  130 CONTINUE
*
*     Write data on disk and print reflections.
*
      JJ=1
      I=0
      NC=1
      DO 160 IS=1,MREFL
      LLS=LL(IS)
      IF(INDEX(IS).EQ.-1) THEN
         I=I+1
         TL(IS)=I
         ICODE(JJ)=JCODE(LLS)
         DO 140 J=1,3
  140    HKL(J,JJ)=IH(J,LLS)
         SJJ(JJ)=M(LLS)*PI*PLOR(LLS)
         JJ=JJ+1
         IOUT(NC)=I
      ELSE
         TL(IS)=0
         IOUT(NC)=8888
      END IF
      ICOUT(NC)=JCODE(LLS)
      DO 150 J=1,3
  150 IHOUT(J,NC)=IH(J,LLS)
      MOUT(NC)=M(LLS)
      HWOUT(NC)=HW(IS)
      TOUT(NC)=T(IS)
      IF(NC.NE.3.AND.IS.NE.MREFL) GO TO 160
      WRITE(7,1012) (IOUT(J),ICOUT(J),(IHOUT(IZ,J),IZ=1,3),MOUT(J),
     *               HWOUT(J),TOUT(J),J=1,NC)
      NC=0
  160 NC=NC+1
*
*     Correct for background and print profile intensities.
*
      CLOSE(UNIT=7,STATUS='KEEP')
      OPEN(UNIT=8,STATUS='NEW',FORM='FORMATTED')
      WRITE(8,1013) IALFA,IOMEG,ISTAP
      WRITE(8,1014)
      I=IALFA
      JJ=1
      IZ=0
      NZ=0
  170 Y=YOBS(JJ)
      IF(Y.EQ.0) THEN
         NZ=NZ+1
         IZ=1
         YOBS(JJ)=0
         WEIGHT(JJ)=0
         GO TO 220
      ELSE
         DO 180 J=1,NB+1
         J1=J-1
         IF(J1.EQ.0) J1=1
         IF(IFIX(POS(J)).LE.I) GO TO 180
         BACK=BAK(J1)
         IF(J.NE.J1)
     *   BACK=BACK+(I-POS(J1))/(POS(J)-POS(J1))*(BAK(J)-BAK(J1))
         GO TO 190
  180    CONTINUE
      END IF
  190 W=250/Y
      YL=Y-BACK
      IF(YL.EQ.0) YL=1
  200 IF(IZ.NE.1) GO TO 210
      IF(NZ.LE.4) THEN
         WRITE(8,1015) (0,K=1,NZ)
      ELSE
         WRITE(8,1015) 0,0
         WRITE(8,1016) NZ-4
         WRITE(8,1015) 0,0
      END IF
      IF(I.GT.IOMEG) GO TO 230
      IZ=0
      NZ=0
  210 YOBS(JJ)=YL
      WEIGHT(JJ)=W
*
*     Write profile on disk file.
*
      K1=MOD(IORD(JJ),OFFSET)
      K2=IORD(JJ)/OFFSET
      K1=TL(K1)+0.5
      K2=TL(K2)+0.5
      IORD(JJ) = K1 + OFFSET * K2
      WRITE(8,1017) I,Y,BACK,YL,W,K1,K2,(STAR,K=K1,K2)
  220 I=I+ISTAP
      JJ=JJ+1
      IF(I.LE.IOMEG) GO TO 170
      IF(IZ.EQ.1) GO TO 200
  230 CLOSE(UNIT=8,STATUS='KEEP')
      NCHAN=(IOMEG-IALFA+ISTAP)/ISTAP
      RETURN

 1000 FORMAT('1PREPARATION OF NEUTRON AND X-RAY POWDER PROFILE INTEN',
     *       'SITIES FOR MAGNETIC AND NUCLEAR STRUCTURE REFINEMENT ',
     *       'PROGRAM.',/,'0FORTRAN 77 VERSION  8-OCT-1987',/)
 1001 FORMAT('0WAVELENGTH =',F7.4)
 1002 FORMAT('0HALFWIDTH PARAMETERS')
 1003 FORMAT('   T =',F10.3,'   U =',F10.3,'   V =',F10.3,'   W =',F10.3,
     *       '   X =',F10.3,'   Y =',F10.3)
 1004 FORMAT('0ZEROPOINT  =',2F6.2)
 1005 FORMAT('   CTHM   Polarization   CUT-value',/,F8.5,F12.5,F12.2)
 1006 FORMAT('0BACKGROUND',/,' POSITION',4X,'INTENSITY')
 1007 FORMAT(1X,F7.0,6X,F6.0)
 1008 FORMAT('0EXCLUDED REGIONS = FROM     TO',/,(20X,I6,I7))
 1009 FORMAT('0NO EXCLUDED REGIONS.')
 1010 FORMAT(/,3('    NO. CODE  H  K  L  MULT   HW     POSN',2X))
 1011 FORMAT(3(12X,3I3,22X))
 1012 FORMAT(3(1X,2I5,I4,2I3,I5,F7.1,F9.1,1X))
 1013 FORMAT('0DIAGRAM FROM',I6,' TO',I6,' IN STEPS OF',I3,
     *       ' 1/100ths DEGREE')
 1014 FORMAT('0  POSN   I+B',6X,'B     I',7X,'W',8X,'PEAKS',7X,'SHAPE')
 1015 FORMAT(1X,I12)
 1016 FORMAT(6X,'(',I4,' ZEROS)')
 1017 FORMAT(1X,I6,3F7.0,F9.4,I5,' -',I4,5X,(79A1))
      END
