* 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