SUBROUTINE MRQSHELD(TOBS,V1OBS,V2V1OBS,SIG1,SIG21,NVPTS,A,CHISQ, 1 STARNAME,DO_1,DO_2,PUNCERT3D,BADPTS) C This subroutine estimates selected parameters for a night's binary C star data using the Levenberg-Marquardt method. C 4 Mar 1991 J.T. Armstrong C 5 Sep 1991 J.T. Armstrong adapted for radial velocities C 23 June 1994 J.T. Armstrong adapted for radial velocity C differences IMPLICIT UNDEFINED (A-Z) INCLUDE 'VELPARM.INC' REAL*8 TOBS(MAXVDATA),V1OBS(MAXVDATA),V2V1OBS(MAXVDATA) REAL*8 SIG1(MAXVDATA),SIG21(MAXVDATA) REAL*8 A(MMAX),COVAR(MMAX,MMAX),ALPHA(MMAX,MMAX) REAL*8 ALAMDA, CHISQ, OCHISQ, DELCHISQ, NCHISQ REAL*8 TCHISQ,RFREE REAL*8 T,A0,A1,B1,B2,XP,EXPCHISQ,DCHISQ REAL*8 DUMMY(MMAX),COVARI,COVARJ,CONF REAL*8 PUNCERT3D(MMAX,-1:1) INTEGER*4 ITST, IERR, ITER, NVPTS, MA, IPARM, I, J INTEGER*4 LISTA(MMAX),MFIT,FREE,BADPTS CHARACTER*1 ANSWER CHARACTER*8 PARMNAME(MMAX) CHARACTER*12 COVFILE CHARACTER*80 STARNAME LOGICAL*1 REPEAT,DO_1,DO_2 DATA A0,A1,B1,B2/2.30753,0.27061,0.99229,0.04481/ DATA PARMNAME/'a (mas) ','eccen ','period ','epoch ', 1 'inclin ','asc node','arg peri','K1 ', 2 'K2 ','v(0) ','dv(0)/dt'/ C Get the list of parameters to vary MFIT = 0 MA = MMAX DO I=1,MMAX LISTA(I) = 0 END DO 1 WRITE(OUTC,*) ' List parameters to vary (2,3,4,7,8,9,10,11):' WRITE(OUTC,*) ' (0 = end, -1 for list)' 2 READ(INC,*) IPARM IF(IPARM.EQ.2 .OR. IPARM.EQ.3 .OR. IPARM.EQ.4 .OR. IPARM.EQ.7 .OR. * IPARM.EQ.8 .OR. IPARM.EQ.9 .OR. IPARM.EQ.10 .OR. IPARM.EQ.11) * THEN MFIT = MFIT+1 LISTA(MFIT) = IPARM GO TO 2 ELSE IF(IPARM.EQ.-1) THEN DO J=1,11 WRITE(OUTC,'(1X,I2,3X,A)') J, PARMNAME(J) END DO GO TO 1 ELSE IF(IPARM.NE.0) THEN WRITE(OUTC,*) ' This parameter number not allowed to vary.' GO TO 1 END IF WRITE(OUTC,*) ' Fractional change in CHISQ for convergence: ' READ(INC,*) DELCHISQ C Get best fit to MFIT parameters C WRITE(OUTC,'(A,I3)') ' MRQSHELL: MFIT = ',MFIT 10 IF (MFIT.GT.1) THEN ALAMDA = -1. ! Initialization step ITER = 0 ITST = 0 D WRITE(OUTC,*) ' Going to MRQMIND first time' CALL MRQMIND(TOBS,V1OBS,V2V1OBS,SIG1,SIG21,NVPTS,A,MA,LISTA, 1 MFIT,COVAR,ALPHA,CHISQ,ALAMDA,DO_1,DO_2) 100 OCHISQ = CHISQ ! Loop to here if no convergence C WRITE(OUTC,*) 'Try again:' ITER = ITER + 1 C WRITE(OUTC,'(A,G11.4)') ' MRQSHELL: OCHISQ = ',OCHISQ WRITE(OUTC,'(A,I5)') ' MRQSHELL: ITER = ',ITER C Call MRQMIND and observe change in CHISQ D WRITE(OUTC,*) ' Going to MRQMIND second time' CALL MRQMIND(TOBS,V1OBS,V2V1OBS,SIG1,SIG21,NVPTS,A,MA,LISTA, 1 MFIT,COVAR,ALPHA,CHISQ,ALAMDA,DO_1,DO_2) D WRITE(OUTC,*) ' Back from MRQMIND second time' IF( .NOT.REPEAT .AND. OCHISQ/CHISQ-1..GT.0. 1 .AND. OCHISQ/CHISQ-1..LT.DELCHISQ ) THEN ITST = ITST + 1 END IF IF (ITST.LT.2 .AND. ITER.LT.20) GO TO 100 IF (ITER.GE.100) THEN WRITE(OUTC,'(A)') ' Iteration limit reached.' END IF C ELSE C CALL CADJMOD C CALL BINVIS C CALL CAGRFAC END IF ALAMDA = 0. ! To get final COVAR and ALPHA D WRITE(OUTC,*) ' Going to MRQMIND last time' CALL MRQMIND(TOBS,V1OBS,V2V1OBS,SIG1,SIG21,NVPTS,A,MA,LISTA, 1 MFIT,COVAR,ALPHA,CHISQ,ALAMDA,DO_1,DO_2) C Write the results to the screen WRITE(OUTC,*) ' Best-fit values:' WRITE(OUTC,'(A,F7.4)') ' ellipticity = ',A(2) WRITE(OUTC,'(A,F9.5)') ' period (days) = ',A(3) WRITE(OUTC,'(A,F11.2)') ' epoch = JD ',A(4) WRITE(OUTC,'(A,F7.2)') ' arg peri (deg) = ',A(7)/RPDEG WRITE(OUTC,'(A,F8.3)') ' K1 (km/s) = ',A(8) WRITE(OUTC,'(A,F8.3)') ' K2 (km/s) = ',A(9) WRITE(OUTC,'(A,F6.2)') ' velocity (km/s) = ',A(10) WRITE(OUTC,'(A,F7.4)') ' accel (km/s/yr) = ',A(11) C Get the confidence level desired WRITE(OUTC,*) ' Desired confidence level (> 50 percent): ' 200 READ(INC,*) CONF IF (CONF.GT.1.) CONF = CONF/100. ! if CONF given in percent IF (CONF.LT.0.5) THEN WRITE(OUTC,*) ' Confidence level must be > 50 percent: ' GO TO 200 END IF C Calculate CHISQ corresponding to indicated confidence level. C See Abramowitz & Stegun 26.2.22 and 26.4.17. T = DSQRT(-2.*LOG(1.-CONF)) XP = T - (A0+A1*T)/(1.+B1*T+B2*T*T) FREE = 2*NVPTS-MFIT-BADPTS RFREE = FREE ! RFREE is REAL*8 NCHISQ = (1.-(2./(9.*RFREE))+XP*DSQRT(2./(9.*RFREE)))**3 C Calculate CHISQ expected for this number of degrees of freedom. C (Same calculation as above, with CONF = 0.5 => XP = 0.) EXPCHISQ = (1.-(2./(9.*RFREE)))**3 DCHISQ = NCHISQ-EXPCHISQ WRITE(OUTC,'(A,I3,A,F5.3,A,F5.3,A,F4.1,A)') 1 ' Expect Chi^2 (',FREE,' d.o.f.) = ',EXPCHISQ,' (min) or ', 2 NCHISQ,' (',CONF*100.,'% confidence).' C Expect Chi^2 (123 d.o.f.) = 0.988 (min) or 1.234 (99.9% confidence). WRITE(OUTC,'(A,F6.3,A,F6.3,A,F4.1,A)') 1 ' Observe Chi^2(min) = ',CHISQ/RFREE,'; seek Chi^2 = ',NCHISQ, 2 ' at ',CONF*100.,'% confidence.' C Observe Chi^2(min) = 1.234; seek Chi^2 = 2.345 at 99.9% confidence. IF(CHISQ/RFREE.GT.EXPCHISQ) THEN NCHISQ = NCHISQ*CHISQ/(EXPCHISQ*RFREE) WRITE(OUTC,'(A,G11.4)') ' Scaled NCHISQ = ',NCHISQ END IF C Calculate the confidence intervals from the diagonal elements of the C covariance matrix (see, e.g., Numerical Recipes, p. 536, eq. 14.5.4). C In the approximation of chi^2 forming a parabolic hypersurface in C parameter space, this will give the uncertainty in each parameter C when all the others are free to adjust to the perturbation. DO I=1,MFIT PUNCERT3D(LISTA(I),0) = 1 SQRT(RFREE*DCHISQ*COVAR(LISTA(I),LISTA(I))) IF(LISTA(I).EQ.5 .OR. LISTA(I).EQ.6 .OR. LISTA(I).EQ.7) 1 THEN WRITE(OUTC,'(A,A,A,G11.4)') 1 ' Uncertainty in ',PARMNAME(LISTA(I)),' = ', 2 PUNCERT3D(LISTA(I),0)/RPDEG ELSE WRITE(OUTC,'(A,A,A,G11.4)') 1 ' Uncertainty in ',PARMNAME(LISTA(I)),' = ', 2 PUNCERT3D(LISTA(I),0) END IF END DO C Repeat with smaller DELCHISQ if desired WRITE(OUTC,'(A,F7.3,A)') ' Convergence is slow (',CHISQ,')' WRITE(OUTC,'(A)') ' Continue (Y or N)?' READ(INC,'(A1)') ANSWER CALL CAPS(ANSWER) IF(ANSWER.EQ.'Y') THEN DELCHISQ = DELCHISQ/5. GO TO 10 END IF C Write the covariance matrix WRITE(OUTC,*) 1 'Covariance matrix to screen (S), file (F), or both (B)? ' READ(INC,'(A1)') ANSWER CALL CAPS(ANSWER) IF(ANSWER.EQ.'S' .OR. ANSWER.EQ.'B') THEN WRITE(OUTC,1000) STARNAME(1:8) WRITE(OUTC,1001) (PARMNAME(LISTA(I)),I=1,MFIT) DO I=1,MFIT WRITE(OUTC,1002) PARMNAME(LISTA(I)), 1 (COVAR(LISTA(I),LISTA(J)),J=1,MFIT) END DO WRITE(OUTC,1003) WRITE(OUTC,1001) (PARMNAME(LISTA(I)),I=1,MFIT) DO I=1,MFIT COVARI=SQRT(COVAR(LISTA(I),LISTA(I))) DO J=1,MFIT COVARJ=SQRT(COVAR(LISTA(J),LISTA(J))) DUMMY(J)=COVAR(LISTA(I),LISTA(J))/(COVARI*COVARJ) END DO WRITE(OUTC,1002) PARMNAME(LISTA(I)), 1 (DUMMY(J),J=1,MFIT) END DO DO I=1,MFIT IF(LISTA(I).EQ.5 .OR. LISTA(I).EQ.6 .OR. LISTA(I).EQ.7) 1 THEN WRITE(OUTC,'(A,A,A,G11.4)') 1 ' Uncertainty in ',PARMNAME(LISTA(I)),' = ', 2 PUNCERT3D(LISTA(I),0)/RPDEG ELSE WRITE(OUTC,'(A,A,A,G11.4)') 1 ' Uncertainty in ',PARMNAME(LISTA(I)),' = ', 2 PUNCERT3D(LISTA(I),0) END IF END DO END IF IF(ANSWER.EQ.'F' .OR. ANSWER.EQ.'B') THEN 300 WRITE(OUTC,*) 'Enter file for covariance matrix output: ' READ(INC,'(A12)') COVFILE OPEN(UNIT=OUTVFILE,FILE=COVFILE,STATUS='NEW',IOSTAT=IERR) IF(IERR.NE.0) THEN WRITE(OUTC,*) 'Cannot open file ',COVFILE GO TO 300 END IF WRITE(OUTVFILE,1000) STARNAME(1:8) WRITE(OUTVFILE,1010) WRITE(OUTVFILE,1011) (PARMNAME(LISTA(I)),I=1,MFIT) WRITE(OUTVFILE,1010) DO I=1,MFIT WRITE(OUTVFILE,1012) PARMNAME(LISTA(I)), 1 (COVAR(LISTA(I),LISTA(J)),J=1,MFIT) WRITE(OUTVFILE,1010) END DO WRITE(OUTVFILE,1010) WRITE(OUTVFILE,1010) WRITE(OUTVFILE,1003) WRITE(OUTVFILE,1010) WRITE(OUTVFILE,1011) (PARMNAME(LISTA(I)),I=1,MFIT) WRITE(OUTVFILE,1010) DO I=1,MFIT COVARI=SQRT(COVAR(LISTA(I),LISTA(I))) DO J=1,MFIT COVARJ=SQRT(COVAR(LISTA(J),LISTA(J))) DUMMY(J)=COVAR(LISTA(I),LISTA(J))/(COVARI*COVARJ) END DO WRITE(OUTVFILE,1012) PARMNAME(LISTA(I)), 1 (DUMMY(J),J=1,MFIT) WRITE(OUTVFILE,1010) END DO DO I=1,MFIT IF(LISTA(I).EQ.5 .OR. LISTA(I).EQ.6 .OR. LISTA(I).EQ.7) 1 THEN WRITE(OUTVFILE,'(A,A,A,G11.4)') 1 ' Uncertainty in ',PARMNAME(LISTA(I)),' = ', 2 PUNCERT3D(LISTA(I),0)/RPDEG ELSE WRITE(OUTVFILE,'(A,A,A,G11.4)') 1 ' Uncertainty in ',PARMNAME(LISTA(I)),' = ', 2 PUNCERT3D(LISTA(I),0) END IF END DO CLOSE(OUTVFILE) END IF IF(ANSWER.EQ.'S') THEN WRITE(OUTC,*) 'Type to continue.' READ(INC,'(A1)') ANSWER END IF RETURN 1000 FORMAT(' Covariance matrix for ',A8) 1001 FORMAT(1X,5X, 20(4X,A3,4X)) 1002 FORMAT(1X,A3,2X,20(G9.3,2X) ) 1003 FORMAT(' Correlation matrix') 1010 FORMAT(' ') 1011 FORMAT(1X,5X, 20(4X,A3,5X)) 1012 FORMAT(1X,A3,2X,20(G10.4,2X) ) END