*     ***************
      SUBROUTINE ITER
*     ***************
*
      INTEGER OFFSET
      PARAMETER (OFFSET = 2**16)

      LOGICAL VERT
      CHARACTER ATEXT*4,TEXT*80
      COMPLEX CWERF,VALUE,Z
      INTEGER I1(3),ICALC(10),IOBS(10)
      REAL DUMP(999),SOMEGA(999),DERIV(51)
      INCLUDE 'COMMONS.FOR'
      COMMON/CHARS/ATEXT(30),TEXT
      COMMON/DERMAT/DERSTO(999,51)
      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/IO/EPS,RELAX,RELAXB,RELAXC,SUMOBS,SUMW,SUMWBS,SW,IPONS,
     *          LMAT,LOUT,NUM
      SAVE
      DATA SQRPI/1.7724538509/,CONST/11459.155903/,SKEW/0.0/
      DATA SQRLN2/0.832554611158/,GAUSS/0.9394372787/
*
      WRITE(6,1000)
      WRITE(6,1001) ICYCLE
      IF(FORCED) WRITE(6,1002)
      SUMW=0
      SUMCAL=0
      SUMDIF=0
      SUMFOB=0
      SUMDFF=0
      FNUCDI=0
      FNUCOB=0
      FMAGDI=0
      FMAGOB=0
      DO 10 I=1,MAXS+MAXCON
      V(I)=0
      DO 10 J=1,MAXS+MAXCON
   10 SMM(I,J)=0
      YES=(IPONS.EQ.2.OR.IPONS.EQ.3).AND.ICYCLE.EQ.NCYCLE
      TYPE=(LOUT.EQ.1).AND.(ICYCLE.EQ.NCYCLE)
      PUNCH=((IPONS.EQ.1).OR.(IPONS.EQ.3)).AND.(ICYCLE.EQ.NCYCLE)
      DO 20 I=1,NREFL
      IF(.NOT.STEADY) CALL COORD(HKL(1,I),I)
      SOMEGA(I)=0
      FESD(I)=0
   20 FOBS(I)=0
      IF(.NOT.STEADY) CALL PROFIL
      IF(YES) THEN
         OPEN(UNIT=13,STATUS='NEW',FORM='FORMATTED'
     *        ,CARRIAGECONTROL='LIST'                         ! VAX/VMS
     *        )
         WRITE(13,1003) TEXT
         WRITE(13,1004) 0.01*(IALFA-PINSTR(9)),0.01*ISTAP,NCHAN
      END IF
      IF(TYPE) THEN
         IOMEG=IALFA+(NCHAN-1)*ISTAP
         WRITE(6,1005) IALFA,IOMEG,ISTAP
      END IF
      CALL CALCUL
      NO=0
      DO 110 M=1,NCHAN
      ISS=IALFA+(M-1)*ISTAP
      YCALC=0
      IF(ABS(YOBS(M)).LT.0.001) GO TO 100
      DO 40 J=1,MAXS
   40 DERIV(J)=0
*
*     Calculate the contribution of the reflections ORD1 to ORD2 to the
*     derivatives w.r.t. the profile intensity YOBS.
*
      IORD1=MOD(IORD(M),OFFSET)
      IORD2=IORD(M)/OFFSET
      DO 70 J=IORD1,IORD2
      VERT=BRAGG(J).LE.LIM
      TANTH=TAN(BRAGG(J)/CONST)
      DELTA=ISS-BRAGG(J)-PINSTR(9)-PINSTR(10)/TANTH
      DELT=DELTA*DELTA
      CORR=1
      COS2TH=1/(1+TANTH**2)
      SINTH=TANTH*SQRT(COS2TH)
      IF(VERT) THEN
         ASYMM=DELT*SIGN(1.,DELTA)
         CORR=1-PINSTR(1)*ASYMM/TANTH
         IF(CORR.LE.0) CORR=0.0001
      END IF
*
*     Explanation:
*        GHALF  = Gaussian halfwidth 2w(g)
*        BETAG  = Gaussian integral breadth = GHALF/0.9394372787...
*        Y      = Complex part of Voigt function argument
*        DLNDTH = Derivative of LOG(BETAG) with respect to Bragg angle
*        DYDTH  = Derivative of y with respect to Bragg angle
*
*     Function EXPERF(Y) evaluates EXP(Y**2) * ERFC(Y).
*
      GHALF=SQRT(PINSTR(3)/TANTH**2+TANTH*(PINSTR(4)*TANTH+PINSTR(5))+
     *           PINSTR(6))
      Y=SQRLN2*(PINSTR(7)*TANTH+PINSTR(8)/SQRT(COS2TH))/GHALF
      DLNDTH=(-PINSTR(3)/TANTH**3+PINSTR(4)*TANTH+PINSTR(5)/2)
     *       /(COS2TH*GHALF**2)
      DYDTH=SQRLN2*(PINSTR(7)+PINSTR(8)*SINTH)/(GHALF*COS2TH)-Y*DLNDTH
      X=SQRPI*GAUSS*DELTA/GHALF
      Z=CMPLX(X,Y)
*
*     Voigtian profile shape function - Gautschi algorithm CWERF.
*
      VALUE=SJJ(J)*CORR*GAUSS*CWERF(Z)/GHALF
      OMEGA=REAL(VALUE)
      SOMEGA(J)=SOMEGA(J)+OMEGA
      DUMP(J)=OMEGA*(FNUC(J)+FMAG(J))
      YCALC=YCALC+DUMP(J)
      TERM1=X*OMEGA-Y*AIMAG(VALUE)
      TERM2=X*AIMAG(VALUE)+Y*OMEGA-SJJ(J)*CORR*GAUSS/(SQRPI*GHALF)
*
*     DER = OMEGA for those parameters that are not part of the
*     argument of the profile shape function.
*
      DO 70 K=1,MAXS
      DER=OMEGA
*
*     Cell constant derivatives. Note that CONST includes value 2
*     and thus CONST may not be multiplied with it any more.
*     Note also that zeropoint angle dependence affects this term.
*
      DO 60 L=1,6
      IF(LCRYST(L).EQ.K) THEN
         DER=2*TERM1*(CONST*SQRPI*GAUSS*(1+PINSTR(10)*(1+TANTH**2)/2)
     *       /GHALF-X*DLNDTH)+2*TERM2*DYDTH-OMEGA*DLNDTH
         GO TO 70
      END IF
   60 CONTINUE
*
*     Gaussian halfwidth parameter derivatives.
*
      IF(LINSTR(3).EQ.K.OR.LINSTR(4).EQ.K.OR.LINSTR(5).EQ.K.OR.
     *   LINSTR(6).EQ.K) THEN
         DER=(2*X*X-2*Y*Y-1)*OMEGA-4*X*Y*AIMAG(VALUE)+
     *        2*Y*SJJ(J)*CORR*GAUSS/(SQRPI*GHALF)
*
*     Lorentzian halfwidth parameter derivatives.
*
      ELSE IF(LINSTR(7).EQ.K.OR.LINSTR(8).EQ.K) THEN
         DER=2*TERM2
*
*     Zeropoint derivative - constant part.
*
      ELSE IF(LINSTR(9).EQ.K) THEN
         IF(VERT) THEN
            HALF=GHALF*APPROX(Y)/EXPERF(Y)
            X1=PINSTR(1)*SIGN(1.0,DELTA)*HALF**2/CORR
         ELSE
            X1=0
         END IF
         DER=2*(1+X1)*TERM1
*
*     Zeropoint derivative - cotangent part.
*
      ELSE IF(LINSTR(10).EQ.K) THEN
         IF(VERT) THEN
            HALF=GHALF*APPROX(Y)/EXPERF(Y)
            X1=PINSTR(1)*SIGN(1.0,DELTA)*HALF**2/CORR
         ELSE
            X1=0
         END IF
         DER=2*(1+X1)*TERM1/TANTH
*
*                                  Asymmetry parameter derivative.
*
      ELSE IF(LINSTR(1).EQ.K.AND.VERT) THEN
         DER=ASYMM*OMEGA/CORR
      END IF
*
   70 DERIV(K)=DERSTO(J,K)*DER+DERIV(K)
*
*     Form sums.
*
      YCALC=YCALC/SCALE
      DELTA=YOBS(M)-YCALC
      IF(YCALC.NE.0) THEN
         DO 80 J=IORD1,IORD2
         X=DUMP(J)/YCALC
         IF(TYPE) FESD(J)=X*X/WEIGHT(M)+FESD(J)
   80    FOBS(J)=FOBS(J)+X*YOBS(M)
         SUMCAL=SUMCAL+YCALC
      END IF
      SUMW=SUMW+WEIGHT(M)*DELTA*DELTA
      SUMDIF=SUMDIF+ABS(DELTA)
      DELTA=DELTA*SCALE
*
*     Sum normal matrix elements and vector.
*
      DO 90 J=1,MAXS
      X=WEIGHT(M)*DERIV(J)
      V(J)=V(J)+DELTA*X
      DO 90 K=J,MAXS
   90 SMM(J,K)=SMM(J,K)+X*DERIV(K)
  100 IF(.NOT.(YES.OR.TYPE)) GO TO 110
      NO=NO+1
      ICALC(NO)=YCALC+0.5
      IOBS(NO)=NINT(YOBS(M))
      IF(NO.NE.10.AND.M.LT.NCHAN) GO TO 110
*
*     Print (TYPE=TRUE) and punch (YES=TRUE) calculated and observed
*     profile intensities.
*
      IF(TYPE) WRITE(6,1006) ISS-(NO-1)*ISTAP,(IOBS(J),ICALC(J),J=1,NO)
      IF(YES) WRITE(13,1007) (IOBS(J),J=1,NO)
      IF(YES) WRITE(13,1007) (ICALC(J),J=1,NO)
      NO=0
  110 CONTINUE
      IF(YES) WRITE(13,1008)
      DO 120 I=1,NREFL
      IF(FOBS(I).LE.0) THEN
         FOBS(I)=0
         FESD(I)=0
      END IF
      BB=SJJ(I)/(ISTAP*SCALE)
      IF(TYPE) FESD(I)=BB*SQRT(FESD(I)*250)/SOMEGA(I)
      FNUC(I)=FNUC(I)*BB
      FMAG(I)=FMAG(I)*BB
  120 FOBS(I)=FOBS(I)*BB/SOMEGA(I)
*
*     Update all variables and print the data.
*
      CALL UPDATE(SKEW)
*
*     On last cycle all selected data is printed.
*
      IF(.NOT.TYPE) GO TO 170
      WRITE(6,1009)
      WRITE(6,1011) NH
      DO 160 I=1,NREFL
      K=FNUC(I)+FMAG(I)+0.5
      L=NINT(FOBS(I)-K)
      TANTH=TAN(BRAGG(I)/CONST)
      M=BRAGG(I)+PINSTR(9)+PINSTR(10)/TANTH+0.5
      DO 130 J=1,3
  130 I1(J)=NINT(HKL(J,I)*NH(J))
      INUC=FNUC(I)+0.5
      IMAG=FMAG(I)+0.5
      KOBS=NINT(FOBS(I))
      IESD=FESD(I)+0.5
      GHALF=SQRT(PINSTR(3)/TANTH**2+TANTH*(PINSTR(4)*TANTH+PINSTR(5))+
     *           PINSTR(6))
      CAUCHY=PINSTR(7)*TANTH+PINSTR(8)*SQRT(1+TANTH**2)
      Y=SQRLN2*CAUCHY/GHALF
      TOTAL=GHALF*APPROX(Y)/EXPERF(Y)
      WRITE(6,1011) I1,M,INUC,IMAG,K,KOBS,L,IESD,TEXCOR(I),GHALF,CAUCHY,
     *              TOTAL,Y
  160 CONTINUE
  170 DO 180 I=1,NREFL
      IF(FOBS(I).EQ.0.AND.FNUC(I)+FMAG(I).EQ.0) GO TO 180
      SUMDFF=SUMDFF+ABS(FOBS(I)-ABS(FNUC(I)+FMAG(I)))
      SUMFOB=SUMFOB+FOBS(I)
      X=FOBS(I)/(FNUC(I)+FMAG(I))
      FNUCDI=FNUCDI+ABS(FNUC(I)-FNUC(I)*X)
      FNUCOB=FNUCOB+FNUC(I)*X
      FMAGDI=FMAGDI+ABS(FMAG(I)-FMAG(I)*X)
      FMAGOB=FMAGOB+FMAG(I)*X
  180 CONTINUE
*
*     Calculate and print sums.
*
      WRITE(6,1012) 100*SUMDFF/SUMFOB,100*SUMDIF/SUMOBS,
     *              100*SQRT(SUMW/SUMWBS)
      WRITE(6,1013) 100*SQRT(250*(NUM-MAXS+MAXCON)/SUMWBS)
      WRITE(6,1014) 100*FNUCDI/FNUCOB
      IF(HALPER.OR.UNIAX.OR.CUBIC) WRITE(6,1015) 100*FMAGDI/FMAGOB
      I=NUM-MAXS+MAXCON
      WRITE(6,1016) I
      WRITE(6,1017) SUMDIF,SUMOBS,SUMCAL,SUMWBS/250,SUMDFF,SUMFOB,
     *              FNUCDI,FNUCOB,FMAGDI,FMAGOB,SUMW/(250*I),SKEW
*
      CALL TIMES(1)
      RETURN
 1000 FORMAT('0',128('+'))
 1001 FORMAT(/,'0CYCLE NUMBER',I3)
 1002 FORMAT(' *FORCED TERMINATION*')
 1003 FORMAT(A80)
 1004 FORMAT(2F10.4,I10)
 1005 FORMAT('0DIAGRAM IS FROM',I6,' TO',I6,' IN STEPS OF',I3,
     *       ' *1/100THS DEGREE',//,' POSITION',10(' YOBS YCALC '))
 1006 FORMAT(1X,21I6)
 1007 FORMAT(10I8)
 1008 FORMAT(3X,'-1000',/,'  -10000')
 1009 FORMAT(/,28X,'I N T E N S I T I E S',39X,'H A L F W I D T H S',/,
     *   '   H  K  L  POSITION  NUCLEAR MAGNETIC   TOTAL  OBSERVED DI',
     *   'FFERENCE  ESD  TEXTURE  GAUSSIAN LORENTZIAN   TOTAL',8X,'Y')
 1011 FORMAT(1X,3I3,I8,I10,3I9,I10,I7,F9.5,4F10.4)
 1012 FORMAT('0R-FACTORS=',3F8.2)
 1013 FORMAT(' EXPECTED =',8X,F8.2)
 1014 FORMAT('0R-FACTOR(NUCLEAR) =',F7.2)
 1015 FORMAT(' R-FACTOR(MAGNETIC)=',F7.2)
 1016 FORMAT('0N-P+C=',I8)
 1017 FORMAT(/,7X,'SUMYDIF',7X,'SUMYOBS',6X,'SUMYCALC',6X,'SUMWYOBSSQ',
     1       5X,'SUMIDIF',7X,'SUMIOBS',/,1X,6E14.4,//,6X,'SUMNUCDIF',5X,
     2       'SUMNUCOBS',5X,'SUMMAGDIF',5X,'SUMMAGOBS',6X,'RESIDUAL',4X,
     3       'LG(SKEWNESS)',/,1X,5E14.4,F16.8,/)
      END
