
*     Least squares fit with conditions
 
*
*     CONDFIT programs
*
      SUBROUTINE SIMNIT(NI,NJ,JDEBUG,EPSIL)
*     initialize and define dimension  NI of X, NJ of Y/F, and
*                    define debug flag JDEBUG, and epsilon (F)
*
      PARAMETER (LUNP=6)
*     =========================== MACRO ================================
*     Parameter statement for basic dimensions
      PARAMETER (NXMAX=50,NFMAX=20)
      PARAMETER (NDIMR=NXMAX+NFMAX)
*     Definition of common for APLCON/ERRPRP/SIM... subroutines
      COMMON/SIMCOM/NX,MYF,NUM,IFLG,INIT,LUNSIM,IDEBUG,ND,CHSQ,EPSF,
     +      ISTAT,
     +      XD(2),XL(2,NXMAX),ST(NXMAX),FC(NXMAX),H(NDIMR)
*     =========================end=of=macro=============================
      COMMON/MATCOM/A(1000)
*     ...
      NX    =NI
      MYF   =NJ
      LUNSIM=LUNP
      IDEBUG=JDEBUG
      ND    =0
      EPSF  =EPSIL
      INIT  =0
      ISTAT =0
*     clear derivative matrix A and ...
      DO 10 IJ=1,NX*MYF
   10 A(IJ)=0.0
*     ... limits XL and steps ST
      DO 20 I=1,NX
      XL(1,I)=0.0
      XL(2,I)=0.0
   20 ST(  I)=0.0
  100 RETURN
      END

      SUBROUTINE SIMSTP(IA,STEP)
*     define step size for num.dif. of variable X(IA)
*
*     =========================== MACRO ================================
*     Parameter statement for basic dimensions
      PARAMETER (NXMAX=50,NFMAX=20)
      PARAMETER (NDIMR=NXMAX+NFMAX)
*     Definition of common for APLCON/ERRPRP/SIM... subroutines
      COMMON/SIMCOM/NX,MYF,NUM,IFLG,INIT,LUNSIM,IDEBUG,ND,CHSQ,EPSF,
     +      ISTAT,
     +      XD(2),XL(2,NXMAX),ST(NXMAX),FC(NXMAX),H(NDIMR)
*     =========================end=of=macro=============================
*     ...
      IF(IA.LT.1.OR.IA.GT.NX) GOTO 100
*     step for numerical differentiation of X(IA)
      ST(IA)=ABS(STEP)
  100 RETURN
      END

      SUBROUTINE SIMNCH(NDEG,CHISQ)
*     get (after convergence) number of degrees of freedom NDEG, and
*                             chi square CHISQ
*
*     =========================== MACRO ================================
*     Parameter statement for basic dimensions
      PARAMETER (NXMAX=50,NFMAX=20)
      PARAMETER (NDIMR=NXMAX+NFMAX)
*     Definition of common for APLCON/ERRPRP/SIM... subroutines
      COMMON/SIMCOM/NX,MYF,NUM,IFLG,INIT,LUNSIM,IDEBUG,ND,CHSQ,EPSF,
     +      ISTAT,
     +      XD(2),XL(2,NXMAX),ST(NXMAX),FC(NXMAX),H(NDIMR)
*     =========================end=of=macro=============================
*     ...
      NDEG =ND
      CHISQ=CHSQ
  100 RETURN
      END

      SUBROUTINE SIMLIM(IA,XLOW,XHIG)
*     define limits for variable X(IA)
*
*
*     =========================== MACRO ================================
*     Parameter statement for basic dimensions
      PARAMETER (NXMAX=50,NFMAX=20)
      PARAMETER (NDIMR=NXMAX+NFMAX)
*     Definition of common for APLCON/ERRPRP/SIM... subroutines
      COMMON/SIMCOM/NX,MYF,NUM,IFLG,INIT,LUNSIM,IDEBUG,ND,CHSQ,EPSF,
     +      ISTAT,
     +      XD(2),XL(2,NXMAX),ST(NXMAX),FC(NXMAX),H(NDIMR)
*     =========================end=of=macro=============================
*     ...
      IF(IA.LT.1.OR.IA.GT.NX) GOTO 100
*     lower or upper limit of X(IA)
      XL(1,IA)=AMIN1(XLOW,XHIG)
      XL(2,IA)=AMAX1(XLOW,XHIG)
  100 RETURN
      END

      SUBROUTINE SIMMAT(VX)
*     check derivative matrix, set flag NUM ...
*        and define steps for num. diff. from covariance matrix
*        NUM = 0      analytical derivatives
*        MUN = 1, 2   numerical derivatives
*                     (=2 for comparison analytica/numerical der.)
*
*     =========================== MACRO ================================
*     Parameter statement for basic dimensions
      PARAMETER (NXMAX=50,NFMAX=20)
      PARAMETER (NDIMR=NXMAX+NFMAX)
*     Definition of common for APLCON/ERRPRP/SIM... subroutines
      COMMON/SIMCOM/NX,MYF,NUM,IFLG,INIT,LUNSIM,IDEBUG,ND,CHSQ,EPSF,
     +      ISTAT,
     +      XD(2),XL(2,NXMAX),ST(NXMAX),FC(NXMAX),H(NDIMR)
*     =========================end=of=macro=============================
      COMMON/MATCOM/A(1000)
      REAL VX(*)
*     ...
      NUM=0
      DO 10 IJ=1,NX*MYF
      IF(A(IJ).NE.0.0) GOTO 20
   10 CONTINUE
*     set NUM flag for zero derivative matrix
      NUM=1
*     force numerical derivatives for first call
   20 IF(NUM.EQ.0.AND.IFLG.EQ.0.AND.IDEBUG.GE.0) NUM=2
      IF(NUM.EQ.0) GOTO 100
*     define steps from covariance matrix for numer. differentiation
      II=0
      DO 30 I=1,NX
      II=II+I
      VII=ABS(VX(II))
      IF(VII.NE.0.0) THEN
         IF(ST(I).NE.0.0) THEN
            ST(I)=AMIN1(ST(I),0.5*SQRT(VII))
         ELSE
            ST(I)=0.5*SQRT(VII)
         END IF
      END IF
   30 CONTINUE
  100 RETURN
      END

      SUBROUTINE SIMDER(X,F,JRET)
*     derivative calculation
*
*     =========================== MACRO ================================
*     Parameter statement for basic dimensions
      PARAMETER (NXMAX=50,NFMAX=20)
      PARAMETER (NDIMR=NXMAX+NFMAX)
*     Definition of common for APLCON/ERRPRP/SIM... subroutines
      COMMON/SIMCOM/NX,MYF,NUM,IFLG,INIT,LUNSIM,IDEBUG,ND,CHSQ,EPSF,
     +      ISTAT,
     +      XD(2),XL(2,NXMAX),ST(NXMAX),FC(NXMAX),H(NDIMR)
*     =========================end=of=macro=============================
      COMMON/MATCOM/A(1000)
      REAL X(*),F(*)
      LOGICAL LIMDEF
      DATA I/0/,XSAVE/0/,ILR/0/
*     ...
*     initialize derivative calculation---------------------------------
      IF(INIT.EQ.0) THEN
*        start with first variable
         I=0
         INIT=1
         GOTO 40
      END IF
      JRET=-1
*     derivative calculation -------------------------------------------
      IF(I.LT.0) THEN
*        another step required, save constraint values ...
         DO 10 J=1,MYF
   10    H(J)=F(J)
*        ... and set next step
         I=-I
         X(I)=XD(2)
         GOTO 100
      END IF
      X(I)=XSAVE
      IJ=I
      DO 20 J=1,MYF
*     calculation of numerical derivative
      IF(ILR.EQ.0) THEN
*        symmetric formula
         DER=0.5*(H(J)-F(J))/ST(I)
      ELSE
*        asymmetric formula
         DER=0.5*(3.0*FC(J)+F(J)-4.0*H(J))/ST(I)
         IF(ILR.EQ.2) DER=-DER
      END IF
*     compare derivatives for NUM=2
      IF(NUM.EQ.2) THEN
         IF(ABS(A(IJ)-DER).GT.0.005*(ABS(A(IJ))+ABS(DER))) THEN
C           WRITE(LUNSIM,101) J,I,A(IJ),DER
         END IF
      END IF
*     insert into A
      A(IJ)=DER
   20 IJ=IJ+NX
*     test end condition
   30 IF(I.EQ.NX) THEN
         JRET=0
         INIT=0
         IF(NUM.EQ.2) NUM=0
         GOTO 100
      END IF
*     next variable
   40 JRET=-1
      I=I+1
      IF(ST(I).EQ.0.0) GOTO 30
      IREDUC=0
      XSAVE=X(I)
*     central differentiation
   50 ILR=0
      LIMDEF=XL(1,I).NE.XL(2,I)
      XD(1)=XSAVE+ST(I)
      IF(LIMDEF.AND.XD(1).GT.XL(2,I)) THEN
*        above upper limit
         XD(1)=XSAVE-ST(I)
         IF(LIMDEF.AND.XD(1).LT.XL(1,I)) GOTO 90
         XD(2)=XSAVE-ST(I)-ST(I)
         IF(LIMDEF.AND.XD(2).LT.XL(1,I)) GOTO 90
*        left differentiation
         ILR=1
      ELSE
         XD(2)=XSAVE-ST(I)
         IF(LIMDEF.AND.XD(2).LT.XL(1,I)) THEN
*           below lower limit
            XD(2)=XSAVE+ST(I)+ST(I)
            IF(LIMDEF.AND.XD(2).GT.XL(2,I)) GOTO 90
*           right differentiation
            ILR=2
         END IF
      END IF
*     first step
      X(I)=XD(1)
      I=-I
      GOTO 100
*     reduce step size
   90 IF(IREDUC.GE.4) GOTO 30
      ST(I)=ST(I)/3.0
      IREDUC=IREDUC+1
      GOTO 50
  100 CONTINUE
      RETURN
  101 FORMAT(/' Derivative dF(',I2,')/dX(',I2,') = ',G15.5,' versus ',
     + G15.5,' (=numerical with >1% deviation) '/)
      END

      SUBROUTINE APLCON(X,VX,F,IRET)
*     Apply constraints F(J) = function of X   J=1,NF
*     to data X(1)...X(nx) with covariance matrix VX.
*     VX(.) in symmetric storage mode
*        1,1  1,2  2,2  1,3  2,3  3,3 ...
*
*     Restriction because of limited dimensions of
*     internal arrays:
*        NX = 50  is maximum
*        NF = 20  is maximum
*
*     Usage
*     =====
*
*        CALL SIMNIT(NX,NF,JDEBUG,EPSILON)
*
*                    JDEBUG = 0   no printout
*                           > 0   more and more printout
*                    EPSILON=     precision required for constraints
*  
*        Now the variables X(1) ... X(NX) have to be defined
*        and their covariance matrix V(1) ...
*        variables may be measured ones or unmeasured ones.
*        unmeasured variables are characterized by zero elements
*        of the covariance matrix (at least the corresponding
*        diagonal element has to be zero). For unmeasured variables
*        the value of x has to be some reasonable initial value.
*        in addition it is necessary to define some step size
*        for numerical differentiation for the unmeasured
*        variables (for measured values the necessary step size
*        is taken from the standard deviation in the covar.
*        matrix). The call to define step ST for variable X(i) is
*           CALL SIMSTP(I,ST)
*        The user may optionally define limits for physical regions
*        of variable x(i) by
*           CALL SIMLIM(I,XLOW,XHIG)
*        The following coding example represents the loop, which
*        performs the constrained fit. during the loop the
*        variables X(1)...X(NX) are modified, to perform the
*        numerical differentiation and to apply corrections
*        during the fit. The user has to supply the code to
*        calculate the constraint equations F(1)...F(NF).
*        Finally (at convergence) the covariance matrix is modified
*        to the matrix for the fitted values, which (for sufficient
*        constraints) has nonzero elements for the previously
*        unmeasured variables.
*
*     10 F(1)=function of X(1) ... X(NX)
*        ...
*        F(NF) = function of X(1) ... X(NX)
*        CALL APLCON(X,VX,F,IRET)
*        IF(IRET.LT.0) GOTO 10
*        CALL SIMNCH(NDEG,CHISQ)      optional to get chisquare
*
*        IRET = 0   convergence reached
*        IRET = 1   bad convergence         no convergence
*        IRET = 2   too many iterations         - " -
*        IRET = 3   unphysical region           - " -
*        IRET = 4   ND less or equal zero       - " -
*
*        Non-convergence is assumed for
*        chisquare above 4 standard deviations for first 10 iterations
*        chisquare above 3 standard deviations for next 10 iterations
*        more than 20 iterations
*
*        The user may scale up his covariance matrix to force more
*        cases with convergence.
*
*        Remark to precision of the constraints:
*        a necessary condition for convergence is the reduction of
*        the constraints to about EPSILON. 
*        Convergence may be
*        difficult due to roundoff-errors, if the required
*        accuracy is too high.
*
*        Remark to differentiation:
*        By default numerical differentiation is done, to calculate
*        the nx*nf elements of the derivative matrix for the
*        constraits w.r.t the variables. Usually this works well.
*        The user may calculate the elements in his program, for
*        reasons of speed or in cases, where the num. diff. fails.
*        The derivative matrix is array A in
*        COMMON/MATCOM/A(.)
*        with the definition
*        df(j)/dx(i) = A(i+nx*(j-1))
*        and has to be defined by the user before the call of
*        APLCON. In the first iteration of the first case the
*        program will automatically compare the elements with
*        values calculated numerically and will printout elements
*        with disagreement.
*
      REAL X(*),VX(*),F(*)
*     =========================== MACRO ================================
*     Parameter statement for basic dimensions
      PARAMETER (NXMAX=50,NFMAX=20)
      PARAMETER (NDIMR=NXMAX+NFMAX)
*     Definition of common for APLCON/ERRPRP/SIM... subroutines
      COMMON/SIMCOM/NX,MYF,NUM,IFLG,INIT,LUNSIM,IDEBUG,ND,CHSQ,EPSF,
     +      ISTAT,
     +      XD(2),XL(2,NXMAX),ST(NXMAX),FC(NXMAX),H(NDIMR)
*     =========================end=of=macro=============================
*     NX       = number of parameters
*     MYF      = number of transformed parameters/constraint equations
*     NUM      = flag for numerical differentiation
*     IFLG     = flag for first case (check derivative)
*     LUNSIM   = printout unit (see parameter statement)
*     ST(50)   = step sizes for numerical differentiation
*     XL(2,50) = lower and upper values of parameters
*     XD(2)    = for current parameter displaced values for differ.
*     H(70)    =
*     FC(50)   = central values of parameters
*     IDEBUG   = debug flag value

*     A        = derivative matrix a/flags during matrix inversion/pulls
*     debug flag = 0   no printout
*                = 1   only warnings are printed (default)
*                = 2,3 more and more printout
*
      PARAMETER (NDIMW=(NDIMR+NDIMR*NDIMR)/2)
      REAL XS(NXMAX),DX(NXMAX),DXP(NXMAX),R(NDIMR),W(NDIMW)
      INTEGER NRD(NXMAX)

*     derivative matrix in MATCOM (see SMINV)
      COMMON/MATCOM/A(1000)
      COMMON/INVCDR/DR(2,100)
      COMMON/CITER/ITER
      CHARACTER*19 TEXT(4)
      DATA         TEXT   /'Chisquare too high ',
     +                     'Too many iterations',
     +                     'Unphysical region  ',
     +                     'ND less or equal 0 '/
      DATA ICNT/0/,          NXF/0/,NCST/0/,FTESTP/0.0/,CHSQP/0.0/
*     statement function for ...
*     chisquare limit for +k sigma and nd degrees of freedom is
*     approximately
      CHLIM(K,ND)=0.5*(FLOAT(K)+SQRT(FLOAT(2*ND+1)))**2
*     ...
      IF(ISTAT.NE.0) GOTO 40
      IFLG=ICNT
      ICNT=ICNT+1
      ISTAT=0
*     initialization----------------------------------------------------
*     define loop parameters
   10 NXF=NX+MYF
      MXF=(NXF*NXF+NXF)/2
*
      CALL SIMMAT(VX)
*     count nr of degrees of freedom
   20 ND=MYF
      II=0
      DO 24 I=1,NX
      II=II+I
*     save initial X values and reset correction DX
      XS(I)=X(I)
      DX(I)=0.0
*     check unphysical region
      IF(XL(1,I).NE.XL(2,I)) THEN
         IF(X(I).LT.XL(1,I).OR.X(I).GT.XL(2,I)) THEN
            IRET=3
            GOTO 99
         END IF
      END IF
      IF(VX(II).LE.0.0) THEN
*        unmeasured variable
         ND=ND-1
*        clear remaining elements
         IJ=II-I
         DO 22 J=1,NX
         IF(J.LE.I) IJ=IJ+1
         VX(IJ)=0.0
         IF(J.GE.I) IJ=IJ+J
   22    CONTINUE
      END IF
   24 CONTINUE
      IF(ND.LE.0) THEN
*        insufficient information
         IRET=4
         GOTO 99
      END IF
      IF(IDEBUG.GE.1) THEN
         WRITE(LUNSIM,101) NX,MYF,ND
         IF(IDEBUG.GE.2) CALL GMPR(X,NX,1,'of initial X-values')
         IF(IDEBUG.GE.3) CALL GMPR(ST,NX,1,'of initial steps')
         IF(IDEBUG.GE.3) CALL SMPRV(X,VX,NX)
      END IF
*     initial value of ftest
      FTESTP=0.0
      DO 26 J=1,MYF
   26 FTESTP=FTESTP+ABS(F(J))
      FTESTP=FTESTP/FLOAT(MYF)
      ITER=0
      NCST=0
      CHSQ=0.0

*     prepare next iteration--------------------------------------------

   30 ISTAT=1
      IRET=-1
*     define right hand side of equation R
      DO 32 I=1,NX
   32 R(I)=0.0
      DO 34 J=1,MYF
*     fc is used in SIMDER
      FC(J)=F(J)
   34 R(NX+J)=-F(J)
      IF(ITER.NE.0) THEN
*        define steps = 0.5 sigma from W
         II=0
         DO 36 I=1,NX
         II=II+I
         IF(ST(I).EQ.0.0.OR.W(II).EQ.0.0) GOTO 36
         ST(I)=AMIN1(ST(I),0.5*SQRT(ABS(W(II))))
   36    CONTINUE
      END IF

*     loop for numerical calculation of derivatives---------------------

   40 IF(ISTAT.NE.1) GOTO 90
      CALL SIMDER(X,F,JRET)
      IF(JRET.LT.0) GOTO 100
      IF(IDEBUG.GE.3.AND.ITER.GE.0)
     +    CALL GMPR(A,MYF,NX,'derivative matrix')

*     construct matrix w and correct vector r---------------------------

*     insert -v and a into w and update r
   50 IJ=(NX*NX+NX)/2
      DO 52 I=1,IJ
   52 W(I)=-VX(I)
      IA=0
      DO 58 J=1,MYF
      DO 54 I=1,NX
      R(NX+J)=R(NX+J)+A(IA+I)*DX(I)
   54 W(IJ+I)=A(IA+I)
      IJ=IJ+NX
      DO 56 K=1,J
      IJ=IJ+1
   56 W(IJ)=0.0
   58 IA=IA+NX

*     calculate step delx-----------------------------------------------

      ITER=ITER+1

*     First part of matrix inversion, making use of
*     the fact, that all elements corresponding to
*     measured variables are already the inverse elements

   60 NM=0
      IF(IDEBUG.GE.4) CALL SMPR(W,NXF,'W before modification ')
      II=0
      DO 61 I=1,NX
      II=II+I
      IF(W(II).LT.0.0) THEN
         DR(1,I)=0.0
         NM=NM+1
         NRD(NM)=I
      ELSE
         W(II)=0.0
         DR(1,I)=1.0
      END IF
   61 CONTINUE
      DO 67 I=NX+1,NXF
      DR(1,I)=1.0
      II=(I*I-I)/2
      DO 63 M=1,NM
      J=NRD(M)
      SUM=0.0
      JK=(J*J-J)/2
      DO 62 K=1,NX
      IF(K.LE.J) JK=JK+1
      IF(DR(1,K).EQ.0.0) SUM=SUM+W(JK)*W(II+K)
   62 IF(K.GE.J) JK=JK+K
   63 H(J)=SUM
      DO 65 K=I,NXF
      IK=(K*K-K)/2
      WJK=0.0
      DO 64 M=1,NM
      J=NRD(M)
   64 WJK=WJK+W(IK+J)*H(J)
   65 W(IK+I)=W(IK+I)+WJK
      DO 66 M=1,NM
      J=NRD(M)
   66 W(II+J)=-H(J)
   67 CONTINUE
*     save right hand side for Chi**2 calculation
      DO 68 J=1,MYF
   68 H(J)=R(NX+J)

*     complete matrix inversion and calculate chisquare

      IF(IDEBUG.GE.4) CALL GMPR(R,1,NXF,'R before solution')
   70 CALL SMINV(W,R,-NXF,1,NRANK)
*     Chi**2 calculation
      CHSQP=CHSQ
      CHSQ=0.0
      DO 72 J=1,MYF
   72 CHSQ=CHSQ-H(J)*R(NX+J)
      IF(IDEBUG.GE.4) CALL SMPR(W,NXF,'W after  SMINV')
      IF(IDEBUG.GE.4) CALL GMPR(R,1,NXF,'R after solution')

      DO 74 I=1,NX
      DXP(I)=DX(I)
   74 DX(I)=R(I)
      GOTO 84
*     make cutstep
   80 DO 82 I=1,NX
   82 DX(I)=0.5*(DX(I)+DXP(I))
*     correct x and return to test constraints
   84 DO 86 I=1,NX
   86 X(I)=XS(I)+DX(I)
      ISTAT=2
*     check unphysical region
      IUNPH=0
      DO 89 I=1,NX
      IF(XL(1,I).NE.XL(2,I)) THEN
C        WRITE(*,*) 'COMPARE ',I,X(I),XL(1,I),XL(2,I)
         IF(X(I).LT.XL(1,I).OR.X(I).GT.XL(2,I)) IUNPH=1
      END IF
   89 CONTINUE
      IF(IUNPH.NE.0) GOTO 95
      IRET=-2
      GOTO 100

*     convergence check-------------------------------------------------

*     calculate value of constraints and compare
   90 FTEST=0.0
      DO 92 I=1,MYF
   92 FTEST=FTEST+ABS(F(I))
      FTEST=FTEST/FLOAT(MYF)
      FTESTP=AMIN1(FTEST,FTESTP)
      IF(IDEBUG.GE.2) WRITE(LUNSIM,102) ITER,NCST,CHSQ,FTEST,FTESTP
      IUNPH=0
      DO I=1,NX
       IF(XL(1,I).NE.XL(2,I)) THEN
          IF(X(I).LT.XL(1,I).OR.X(I).GT.XL(2,I)) IUNPH=1
       END IF
      END DO 
      IF(IDEBUG.GE.3.AND.NCST.EQ.0) THEN
         CALL GMPR(X,NX,1,'of X-values')
         CALL GMPR(F,MYF,1,' of constraint function values')
      END IF
*     divergence/convergence tests
   95 IF(IUNPH.NE.0.OR.FTEST.GT.1.1*FTESTP+EPSF) THEN
*        divergence, make cut steps
         NCST=NCST+1
         IF(NCST.LT.5) GOTO 80
      ELSE IF(NCST.EQ.0) THEN
         IF((ITER.GE.2.OR.CHSQ.LT.CHLIM(1,ND)).AND.
     1   (FTEST.LT.EPSF.AND.CHSQ-CHSQP.LT.0.1)) THEN
*           convergence
            IF(IDEBUG.GE.1) WRITE(LUNSIM,104)ITER
            IF(IDEBUG.GE.2) CALL GMPR(X,NX,1,'of final X-values')
*           pulls
            II=0
            DO 94 I=1,NX
            II=II+I
            A(I)=0.0
            IF(VX(II).GT.0.0) THEN
               IF(VX(II)-W(II).GT.0.0) A(I)=DX(I)/SQRT(VX(II)-W(II))
            END IF
   94       CONTINUE
            DO 96 I=1,(NX*NX+NX)/2
   96       VX(I)=W(I)
            IF(IDEBUG.GE.2) CALL GMPR(A,NX,1,'of pulls')
            IF(IDEBUG.GE.3) CALL SMPRV(X,VX,NX)
            IRET =0
            ISTAT=0
            GOTO 100
         END IF
      END IF
*     continue with iteration or stop
      NCST=0
      IF(ITER.GE. 2.AND.CHSQ.GT.CHLIM(4,ND))  IRET= 1
      IF(ITER.GT.10.AND.CHSQ.GT.CHLIM(3,ND))  IRET= 1
      IF(ITER.GT.20)                          IRET= 2
      IF(IRET.LT.0) GOTO 30
   99 IF(IDEBUG.GE.1) WRITE(LUNSIM,105) ITER,IRET,TEXT(IRET)
      ISTAT=0

  100 RETURN
  101 FORMAT(/' APLCON-   NX =',I3,'     MF =',I3,'     ND =',I3,
     +    '    -APLCON')
  102 FORMAT(' ITERATION',I3,'.',I1,'   Chisquare =',
     + G12.5,'   FTEST =',2G12.5)
  104 FORMAT(' Convergence after',I3,' iterations')
  105 FORMAT(' No convergence (',I3,' ITER, IRET =',I2,')  ',A)
      END

      SUBROUTINE ERRPRP(X,VX,Y,VY,IRET)
*     error propagation for Y = function of X
*        to data X,VX,NX       NX+MYF<51
*
*        CALL SIMDIM(NX,NY)
*     10 Y(1)=FUNCTION OF X
*        CALL ERRPRP(X,VX,Y,VY,IRET)
*        IF(IRET.LT.0) GOTO 10
*
*
      REAL X(*),VX(*),Y(*),VY(*)

*     =========================== MACRO ================================
*     Parameter statement for basic dimensions
      PARAMETER (NXMAX=50,NFMAX=20)
      PARAMETER (NDIMR=NXMAX+NFMAX)
*     Definition of common for APLCON/ERRPRP/SIM... subroutines
      COMMON/SIMCOM/NX,MYF,NUM,IFLG,INIT,LUNSIM,IDEBUG,ND,CHSQ,EPSF,
     +      ISTAT,
     +      XD(2),XL(2,NXMAX),ST(NXMAX),FC(NXMAX),H(NDIMR)
*     =========================end=of=macro=============================

*     DERIVATIVE MATRIX
      COMMON/MATCOM/A(1000)

      DATA ICNT/0/
*     ...
      IF(ISTAT.NE.0) GOTO 20
      IF(IDEBUG.GE.1) WRITE(LUNSIM,101) NX,MYF
      IF(IDEBUG.GE.2) CALL GMPR(X,NX,1,'of X-values')
      IF(IDEBUG.GE.3) CALL SMPRV(X,VX,NX)
      IFLG=ICNT
      ICNT=ICNT+1
      ISTAT=1

      CALL SIMMAT(VX)

      IRET=-1
      DO 10 J=1,MYF
*     FC IS USED IN SIMDER
   10 FC(J)=Y(J)

*     loop for numerical calculation of derivatives---------------------

   20 CALL SIMDER(X,Y,JRET)
      IF(JRET.LT.0) GOTO 100
      IF(IDEBUG.GE.3) CALL GMPR(A,MYF,NX,'derivative matrix')
      IRET=0
      ISTAT=0
      CALL SMAVAT(VX,A,VY,NX,MYF)
      DO 30 I=1,MYF
   30 Y(I)=FC(I)
      IF(IDEBUG.GE.2) CALL GMPR(Y,MYF,1,'of transformed parameters')
      IF(IDEBUG.GE.3) CALL SMPRV(Y,VY,MYF)

  100 RETURN
  101 FORMAT(/'ERRPRP-   NX =',I3,'     MY =',I3,'   -ERRPRP')
      END

      SUBROUTINE ERPCTS 
*     example program ERRPRP
*       error propagation 
      REAL X(2),VX(3),Y(2),VY(3)
      DATA  X/10.0,5.0/
      DATA VX/0.5, 0.0, 2.0/
*     ...
      WRITE(*,*) 'Test of error propagation' 
      CALL SMPRV(X,VX,2)
      CALL SMPR(VX,2,'Covariance matrix before transformation')
      CALL SIMNIT(2,2,3,1.E-6)
 10   Y(1)=SQRT(X(1)**2+X(2)**2)
      Y(2)=ATAN2(X(2),X(1))
      CALL ERRPRP(X,VX,Y,VY,IRET)
      IF(IRET.LT.0) GOTO 10   
      CALL SMPRV(Y,VY,2)
      CALL SMPR(VY,2,'Covariance matrix after transformation')
      END 

      SUBROUTINE APLCTS
*     example program for APLCON:
*        fit two photons to pi zero mass
*
      REAL A(3),B(3),VA(6),VB(6),X(6),VX(21),F(1)
*            first photon:  E phi theta
      DATA A/ 1.0 , 0.1 , 0.5 /
*            second photon: E phi theta
      DATA B/ 2.0 , 0.2 , 0.4 /
*            covariance matrices
      DATA VA/0.0001 , 0.0 , 0.00002 , 0.0 , 0.0 , 0.00001 /
      DATA VB/0.0002 , 0.0 , 0.00002 , 0.0 , 0.0 , 0.00001 /
*     ... 
      WRITE(*,*) 'Test of constrained fit'
*     x = array of measured values
      DO 10 I=1,3
      X(I  ) = A(I)
   10 X(I+3) = B(I)
*     vx = joint covariance matrix
      DO 12 I=1,21
   12 VX(I)=0.0
      CALL SMTOS(VA,1,VX,1,3)
      CALL SMTOS(VB,1,VX,4,3)
*     6 parameters and 1 constraint equations
      CALL SIMNIT(6,1,3,1.E-3)
*     add four vectors
   20 EN = X(1) + X(4)
      PX = X(1)*SIN(X(2))*COS(X(3)) + X(4)*SIN(X(5))*COS(X(6))
      PY = X(1)*SIN(X(2))*SIN(X(3)) + X(4)*SIN(X(5))*SIN(X(6))
      PZ = X(1)*COS(X(2))           + X(4)*COS(X(5))
*     mass squared by square of four vector
      FM2  = EN**2 - PX**2 - PY**2 - PZ**2
*     constraint = calculated mass**2 - pi zero mass**2
      F(1) = FM2 - 0.135**2
*     call fit program
      CALL APLCON(X,VX,F,IREP)
      IF(IREP.LT.0) GOTO 20
*     print fitted mass value
      WRITE(6,101) SQRT(FM2)
  101 FORMAT(' Fitted mass value is',F7.4,' GeV')
      END















