SUBROUTINE SETELL C-------------------------------------------------------------- C This is called by ELLIPSE, and gets the initial parameters of C the ellipse interactively. Model parameters needed are: C MAJOR AXIS, ECCENTRICITY, POSITION ANGLE OF MAJOR AXIS, C COORDINATES OF CENTER C The calling program must set the value of STATE: C 'N' => no points set C 'C', 'A', OR 'B' => the center, one end of the major axis, or one C end of the minor axis has been set C 'AC' OR 'BC' => both ends of an axis plus the center C have been set C 'AB' => one end of each axis has been set C 'ABC' => all points have been set C The idea in resetting points is that the major axis points take C precedence. Thus when a new center is found and a major axis point C already exists, the new major axis is found, then a new minor axis C is found, even if the center changed because of a change in the minor C axis. C-------------------------------------------------------------- IMPLICIT UNDEFINED (A-Z) C-------------------------------------------------------------- C FOR PLOTTING INTEGER*2 IPLOT C-------------------------------------------------------------- C FOR KEYIN REAL*8 ENDTAG, NAMES(6), VALUES(6) INTEGER*4 MODE C-------------------------------------------------------------- C FOR COMMANDS 'D', 'I', 'U' BELOW REAL*8 DIST, DISTANCE INTEGER*4 I, II C-------------------------------------------------------------- C FOR COMMAND 'C', CASES 'AC','BC','ABC' REAL*8 CRAOLD, CDECOLD INTEGER*4 J C-------------------------------------------------------------- C FOR COMMANDS 'L' AND 'E' INTEGER*4 LISTA(7),MFIT REAL*8 COVAR(7,7) C-------------------------------------------------------------- CHARACTER*60 CHKSTRING INCLUDE 'PARAMETR.INC' INCLUDE 'ELLIPSE.INC' LOGICAL FIRST, HCOPY DATA ENDTAG /'/ '/ DATA NAMES /'A','B','C','TH','CRA','CDEC'/ DATA FIRST /.TRUE./ DATA SHOWDEL /.TRUE./ DATA HCOPY /.FALSE./ D WRITE(PR,*) 'SETELL: NPTS = ', NPTS IF (FIRST.AND.STATE.EQ.'N') THEN A_APP = 0. ARA(1) = 0. ARA(2) = 0. ADEC(1) = 0. ADEC(2) = 0. A_ANGLE = 0. B_APP = 0. BRA(1) = 0. BRA(2) = 0. BDEC(1) = 0. BDEC(2) = 0. CRA = 0. CDEC = 0. E_APP = 0. STATE = 'N' ELLFIT = .FALSE. TRUFIT = .FALSE. END IF FIRST = .FALSE. C---------------------------------------------------------------------- C Start the plotting. Return here after a hardcopy. C---------------------------------------------------------------------- 40 IF(HCOPY) THEN plot_dev='/lpm' IER = PGBEGIN(0,PLOTFILE//PLOT_DEV,1,1) ELSE IER = PGBEGIN(0,'/EGA',1,1) END IF C CALL PGASK(.FALSE.) CALL PGSLW(LINEWD) CALL PGSCF(NFONT) IPLOT = 1 IF (IER.NE.1) GO TO 900 IF (.NOT.ELLFIT.AND.STATE.EQ.'ABC') CALL FITELL C RETURN HERE AFTER EACH COMMAND EXCEPT 'Q' AND 'I' 50 CONTINUE CALL SHOWDATA ! PLOT DATA AND ELLIPSE IF (HCOPY) THEN CALL PGEND HCOPY = .FALSE. GO TO 40 END IF CALL EGA_RESTORE_DEFAULT CALL LOCATE (0,1) CALL WRITE_STRING('Command:') C RETURN HERE AFTER COMMAND 'I', SOMETIMES AFTER 'F' 60 CALL PGCURSE(XR,YR,COMND) CALL CAPS(COMND) C------------------------------------------------------------------ C BEGIN COMMAND INTERPRETATION C------------------------------------------------------------------ IF (COMND.EQ.'A') THEN ! SET ENDS OF MAJOR AXIS IF (STATE.EQ.'N') THEN ! SET A1 (1ST MAJ.AX.PT.) ARA(1) = XR ADEC(1) = YR D CALL CURPOS(OUTC,XR,YR) STATE = 'A' CALL LOCATE(0,1) ELSE IF (STATE.EQ.'A') THEN ! SET A2, C (CENTER) ARA(2) = XR ADEC(2) = YR D CALL CURPOS(OUTC,XR,YR) CRA = 0.5 * (ARA(2) + ARA(1) ) CDEC = 0.5 * (ADEC(2) + ADEC(1)) STATE = 'AC' ELSE IF (STATE.EQ.'B') THEN ! SET A1 ARA(1) = XR ADEC(1) = YR D CALL CURPOS(OUTC,XR,YR) STATE = 'AB' ELSE IF (STATE.EQ.'C') THEN ! SET A1, A2 ARA(1) = XR ADEC(1) = YR D CALL CURPOS(OUTC,XR,YR) ARA(2) = 2.*CRA - ARA(1) ADEC(2) = 2.*CDEC - ADEC(1) STATE = 'AC' ELSE IF (STATE.EQ.'AB') THEN ! SET A2, C, RESET B1, SET B2 ARA(2) = XR ADEC(2) = YR D CALL CURPOS(OUTC,XR,YR) CRA = 0.5 * (ARA(2) + ARA(1) ) CDEC = 0.5 * (ADEC(2) + ADEC(1)) C CALCULATE POINT ON MAJOR AXIS WHERE CURRENT MINOR C AXIS INTERSECTS: CALL PERPEND(ARA,ADEC,BRA,BDEC,RAINT,DECINT,1) C POSITION OF NEW MINOR AXIS: BRA(1) = BRA(1) + (CRA - RAINT ) BDEC(1) = BDEC(1) + (CDEC - DECINT) BRA(2) = 2.*CRA - BRA(1) BDEC(2) = 2.*CDEC - BDEC(1) STATE = 'ABC' ELSE IF (STATE.EQ.'AC') THEN ! RESET A1 OR A2; RESET C C FIND NEAREST END OF MAJOR AXIS: DIST1 = (XR - ARA(1))**2 + (YR - ADEC(1))**2 DIST2 = (XR - ARA(2))**2 + (YR - ADEC(2))**2 IF (DIST1.LT.DIST2) THEN J = 1 ELSE J = 2 END IF C ARA(J) = XR ADEC(J) = YR D CALL CURPOS(OUTC,XR,YR) CRA = 0.5 * (ARA(2) + ARA(1) ) CDEC = 0.5 * (ADEC(2) + ADEC(1)) STATE = 'AC' ELSE IF (STATE.EQ.'BC') THEN ! SET A1, A2; RESET B1, B2 ARA(1) = XR ADEC(1) = YR D CALL CURPOS(OUTC,XR,YR) ARA(2) = 2.*CRA - ARA(1) ADEC(2) = 2.*CDEC - ADEC(1) C DROP PERPENDICULAR FROM ONE END OF CURRENT MINOR C AXIS ONTO MAJOR AXIS. CALL PERPEND(ARA,ADEC,BRA,BDEC,RAINT,DECINT,1) C POSITION OF NEW MINOR AXIS: BRA(1) = BRA(1) + ( CRA - RAINT ) BDEC(1) = BDEC(1) + ( CDEC - DECINT ) BRA(2) = 2.*CRA - BRA(1) BDEC(2) = 2.*CDEC - BDEC(1) STATE = 'ABC' ELSE IF (STATE.EQ.'ABC') THEN ! RESET A1 OR A2; RESET ! B1, B2, C C FIND NEAREST END OF MAJOR AXIS: DIST1 = (XR - ARA(1))**2 + (YR - ADEC(1))**2 DIST2 = (XR - ARA(2))**2 + (YR - ADEC(2))**2 IF (DIST1.LT.DIST2) THEN J = 1 ELSE J = 2 END IF ARA(J) = XR ADEC(J) = YR D CALL CURPOS(OUTC,XR,YR) A_ANGLE = ATAN2( (ARA(2)-ARA(1)) , (ADEC(2)-ADEC(1)) ) C CALCULATE NEW CENTER POSITION CRA = 0.5 * (ARA(2) + ARA(1) ) CDEC = 0.5 * (ADEC(2) + ADEC(1)) C CALCULATE NEW MINOR AXIS ENDPOINTS BRA(1) = CRA - B_APP * COS( A_ANGLE ) BDEC(1) = CDEC + B_APP * SIN( A_ANGLE ) BRA(2) = 2.*CRA - BRA(1) BDEC(2) = 2.*CDEC - BDEC(1) C STATE = 'ABC' END IF C------------------------------------------------------------------ ELSE IF (COMND.EQ.'B') THEN ! SET ENDS OF MINOR AXIS IF (STATE.EQ.'N') THEN ! SET B1 BRA(1) = XR BDEC(1) = YR D CALL CURPOS(OUTC,XR,YR) STATE = 'B' ELSE IF (STATE.EQ.'A') THEN ! SET B1 BRA(1) = XR BDEC(1) = YR D CALL CURPOS(OUTC,XR,YR) STATE = 'AB' ELSE IF (STATE.EQ.'B') THEN ! SET B2, C BRA(2) = XR BDEC(2) = YR D CALL CURPOS(OUTC,XR,YR) CRA = 0.5 * (BRA(1) + BRA(2) ) CDEC = 0.5 * (BDEC(1) + BDEC(2)) STATE = 'BC' ELSE IF (STATE.EQ.'C') THEN ! SET B1, B2 BRA(1) = XR BDEC(1) = YR D CALL CURPOS(OUTC,XR,YR) BRA(2) = 2.*CRA - BRA(1) BDEC(2) = 2.*CDEC - BDEC(1) STATE = 'BC' ELSE IF (STATE.EQ.'AB') THEN ! SET C, A2 BRA(2) = XR BDEC(2) = YR D CALL CURPOS(OUTC,XR,YR) C FIND CENTER CRA = 0.5 * ( BRA(1) + BRA(2) ) CDEC = 0.5 * ( BDEC(1) + BDEC(2) ) C FIND SECOND END OF MAJOR AXIS ARA(2) = 2.*CRA - ARA(1) ADEC(2) = 2.*CDEC - ADEC(1) C DROP PERPENDICULAR FROM ONE END OF CURRENT MINOR C AXIS ONTO MAJOR AXIS. CALL PERPEND(ARA,ADEC,BRA,BDEC,RAINT,DECINT,1) C POSITION OF NEW MINOR AXIS: BRA(1) = BRA(1) + ( CRA - RAINT ) BDEC(1) = BDEC(1) + ( CDEC - DECINT ) BRA(2) = 2.*CRA - BRA(1) BDEC(2) = 2.*CDEC - BDEC(1) C STATE = 'ABC' ELSE IF (STATE.EQ.'AC') THEN ! SET B1, B2 BRA(1) = XR BDEC(1) = YR D CALL CURPOS(OUTC,XR,YR) C DROP PERPENDICULAR ONTO MAJOR AXIS. D WRITE(OUTC,*) ' Going to PERPEND w/ state AC' CALL PERPEND(ARA,ADEC,BRA,BDEC,RAINT,DECINT,1) D WRITE(OUTC,*) ' Back from PERPEND w/ state AC' C POSITION OF NEW MINOR AXIS: BRA(1) = BRA(1) + ( CRA - RAINT ) BDEC(1) = BDEC(1) + ( CDEC - DECINT ) BRA(2) = 2.*CRA - BRA(1) BDEC(2) = 2.*CDEC - BDEC(1) STATE = 'ABC' ELSE IF (STATE.EQ.'BC') THEN ! RESET B1 OR B2; RESET C C FIND NEAREST END OF MINOR AXIS: DIST1 = (XR - BRA(1))**2 + (YR - BDEC(1))**2 DIST2 = (XR - BRA(2))**2 + (YR - BDEC(2))**2 IF (DIST1.LT.DIST2) THEN J = 1 ELSE J = 2 END IF BRA(J) = XR BDEC(J) = YR D CALL CURPOS(OUTC,XR,YR) CRA = 0.5 * ( BRA(1) + BRA(2) ) CDEC = 0.5 * ( BDEC(1) + BDEC(2) ) STATE = 'BC' ELSE IF (STATE.EQ.'ABC') THEN ! RESET B1, B2 C FIND NEAREST END OF MINOR AXIS: DIST1 = (XR - BRA(1))**2 + (YR - BDEC(1))**2 DIST2 = (XR - BRA(2))**2 + (YR - BDEC(2))**2 IF (DIST1.LT.DIST2) THEN J = 1 ELSE J = 2 END IF BRA(J) = XR BDEC(J) = YR D CALL CURPOS(OUTC,XR,YR) C DROP PERPENDICULAR ONTO MAJOR AXIS. CALL PERPEND(ARA,ADEC,BRA,BDEC,RAINT,DECINT,J) C POSITION OF NEW MINOR AXIS: RADUMMY(1) = BRA(J) !THESE LINES FOR DECDUMMY(1) = BDEC(J) !DEBUGGING RADUMMY(2) = RAINT DECDUMMY(2) = DECINT BRA(J) = BRA(J) + ( CRA - RAINT ) BDEC(J) = BDEC(J) + ( CDEC - DECINT ) BRA(1) = BRA(J) BDEC(1) = BDEC(J) BRA(2) = 2.*CRA - BRA(1) BDEC(2) = 2.*CDEC - BDEC(1) STATE = 'ABC' END IF C------------------------------------------------------------------ ELSE IF (COMND.EQ.'C') THEN ! SET LOCATION OF CENTER IF (STATE.EQ.'N'.OR.STATE.EQ.'C') THEN ! SET C CRA = XR CDEC = YR D CALL CURPOS(OUTC,XR,YR) STATE = 'C' ELSE IF (STATE.EQ.'A') THEN ! SET C, A2 CRA = XR CDEC = YR D CALL CURPOS(OUTC,XR,YR) ARA(2) = 2.*CRA - ARA(1) ADEC(2) = 2.*CDEC - ADEC(1) STATE = 'AC' ELSE IF (STATE.EQ.'B') THEN ! SET C, B2 CRA = XR CDEC = YR D CALL CURPOS(OUTC,XR,YR) BRA(2) = 2.*CRA - BRA(1) BDEC(2) = 2.*CDEC - BDEC(1) STATE = 'BC' ELSE IF (STATE.EQ.'AB') THEN ! SET C, A2, B1, B2 CRA = XR CDEC = YR D CALL CURPOS(OUTC,XR,YR) ARA(2) = 2.*CRA - ARA(1) ADEC(2) = 2.*CDEC - ADEC(1) C DROP PERPENDICULAR ONTO MAJOR AXIS. CALL PERPEND(ARA,ADEC,BRA,BDEC,RAINT,DECINT,1) C POSITION OF NEW MINOR AXIS: BRA(1) = BRA(1) + ( CRA - RAINT ) BDEC(1) = BDEC(1) + ( CDEC - DECINT ) BRA(2) = 2.*CRA - BRA(1) BDEC(2) = 2.*CDEC - BDEC(1) STATE = 'ABC' ELSE IF (STATE.EQ.'AC') THEN ! RESET C, A1, A2 CRAOLD = CRA CDECOLD = CDEC CRA = XR CDEC = YR D CALL CURPOS(OUTC,XR,YR) ARA(1) = ARA(1) + ( CRA - CRAOLD ) ADEC(1) = ADEC(1) + ( CDEC - CDECOLD) ARA(2) = 2.*CRA - ARA(1) ADEC(2) = 2.*CDEC - ADEC(1) STATE = 'AC' ELSE IF (STATE.EQ.'BC') THEN ! RESET C, B1, B2 CRAOLD = CRA CDECOLD = CDEC CRA = XR CDEC = YR D CALL CURPOS(OUTC,XR,YR) BRA(1) = BRA(1) + ( CRA - CRAOLD ) BDEC(1) = BDEC(1) + ( CDEC - CDECOLD) BRA(2) = 2.*CRA - BRA(1) BDEC(2) = 2.*CDEC - BDEC(1) STATE = 'BC' ELSE IF (STATE.EQ.'ABC') THEN ! MOVE ENTIRE ELLIPSE CRAOLD = CRA CDECOLD = CDEC CRA = XR CDEC = YR D CALL CURPOS(OUTC,XR,YR) ARA(1) = ARA(1) + ( CRA - CRAOLD ) ADEC(1) = ADEC(1) + ( CDEC - CDECOLD) ARA(2) = 2.*CRA - ARA(1) ADEC(2) = 2.*CDEC - ADEC(1) BRA(1) = BRA(1) + ( CRA - CRAOLD ) BDEC(1) = BDEC(1) + ( CDEC - CDECOLD) BRA(2) = 2.*CRA - BRA(1) BDEC(2) = 2.*CDEC - BDEC(1) STATE = 'ABC' END IF C------------------------------------------------------------------ ELSE IF (COMND.EQ.'D') THEN ! DELETE A POINT FROM FIT C FIND NEAREST POINT DISTANCE = 1.E9 DO I = 1, NPTS DIST = (RA(I) - XR)**2 + (DEC(I) - YR)**2 IF (DIST.LT.DISTANCE) THEN DISTANCE = DIST II = I END IF END DO DELETE(II) = .TRUE. DEL_COUNT = DEL_COUNT + 1 C------------------------------------------------------------------ ELSE IF (COMND.EQ.'E') THEN ! ERROR EST. FOR TRUE ORB. IF (TRUFIT) THEN CALL ERREST(ABEST,LISTA,MFIT,COVAR) ELSE CALL LOCATE (0,2) CALL WRITE_STRING('Fit true orbit first') GO TO 60 END IF C------------------------------------------------------------------ ELSE IF (COMND.EQ.'F') THEN ! FIT APPARENT ELLIPSE IF (STATE.EQ.'ABC') THEN CALL FITELL ELSE CALL LOCATE (0,2) CALL WRITE_STRING('Specify all parameters first') GO TO 60 ! RETURN TO PGCURSE LINE END IF C------------------------------------------------------------------ ELSE IF (COMND.EQ.'G') THEN CALL GRIDCHI C------------------------------------------------------------------ ELSE IF (COMND.EQ.'H') THEN ! MAKE HARDCOPY NPLOTS = NPLOTS + 1 HCOPY = .TRUE. WRITE(PLOTFILE(6:8),'(I3.3)') NPLOTS WRITE(6,*) ' CREATING FILE ', PLOTFILE CALL PGEND GO TO 40 C------------------------------------------------------------------ ELSE IF (COMND.EQ.'I') THEN ! IDENTIFY A DATA POINT DISTANCE = 1.E9 DO I = 1, NPTS DIST = (RA(I) - XR)**2 + (DEC(I) - YR)**2 IF (DIST.LT.DISTANCE) THEN DISTANCE = DIST II = I END IF END DO XR = RA(II) YR = DEC(II) CALL LOCATE (0,2) CALL WRITE_STRING('Observation date '//TSTRING(II)) GO TO 60 ! RETURN TO PGCURSE LINE C------------------------------------------------------------------ ELSE IF (COMND.EQ.'K') THEN ! FIND CHI SQUARED FOR A ! GIVEN SET OF ELEMENTS CALL GETCHI CALL PGASK(.TRUE.) CALL PGADVANCE CALL PGASK(.FALSE.) GO TO 50 C------------------------------------------------------------------ ELSE IF (COMND.EQ.'L') THEN ! LEAST-SQUARES ORBIT FIT IF (STATE.EQ.'ABC') THEN CALL ORBFIT(LISTA,MFIT,COVAR) ELSE IF (VALS(49).NE.-1000..AND. 1 VALS(50).NE.-1000..AND. 2 VALS(51).NE.-1000..AND. 3 VALS(52).NE.-1000..AND. 4 VALS(53).NE.-1000..AND. 5 VALS(54).NE.-1000.) THEN CALL ORBFIT(LISTA,MFIT,COVAR) ELSE CALL LOCATE (0,2) CALL WRITE_STRING 1 ('Specify apparent orbit or true orbital elements first') END IF C------------------------------------------------------------------ ELSE IF (COMND.EQ.'Q') THEN ! QUIT TO PROGRAM LEVEL GO TO 80 C------------------------------------------------------------------ ELSE IF (COMND.EQ.'R') THEN ! ROTATE ELLIPSE A_ANGLE = ATAN2( (XR - CRA) , (YR - CDEC) ) D CALL CURPOS(OUTC,XR,YR) IF (STATE.EQ.'N'.OR.STATE.EQ.'A'.OR.STATE.EQ.'B' 1 .OR.STATE.EQ.'C') THEN ! NO ACTION CRA = CRA ELSE IF (STATE.EQ.'AC') THEN ! ROTATE MAJ AX TO PASS ARA(1) = CRA + A_APP * SIN(A_ANGLE) ! THRU POINT ADEC(1) = CDEC + A_APP * COS(A_ANGLE) ARA(2) = 2.*CRA - ARA(1) ADEC(2) = 2.*CDEC - ADEC(1) ELSE IF (STATE.EQ.'BC') THEN ! ROTATE MIN AX TO PASS A_ANGLE = A_ANGLE - 90.*RPDEG ! THRU POINT BRA(1) = CRA + B_APP * COS(A_ANGLE) BDEC(1) = CDEC - B_APP * SIN(A_ANGLE) BRA(2) = 2.*CRA - BRA(1) BDEC(2) = 2.*CDEC - BDEC(1) ELSE IF (STATE.EQ.'ABC') THEN ! ROTATE MAJ AX TO PASS ARA(1) = CRA + A_APP * SIN(A_ANGLE) ! THRU POINT ADEC(1) = CDEC + A_APP * COS(A_ANGLE) ARA(2) = 2.*CRA - ARA(1) ADEC(2) = 2.*CDEC - ADEC(1) BRA(1) = CRA + B_APP * COS(A_ANGLE) BDEC(1) = CDEC - B_APP * SIN(A_ANGLE) BRA(2) = 2.*CRA - BRA(1) BDEC(2) = 2.*CDEC - BDEC(1) END IF C------------------------------------------------------------------ ELSE IF (COMND.EQ.'S') THEN ! SHOW/UNSHOW DELETED POINTS IF (SHOWDEL.EQ..TRUE.) THEN SHOWDEL = .FALSE. ELSE SHOWDEL = .TRUE. END IF C------------------------------------------------------------------ ELSE IF (COMND.EQ.'T') THEN ! FIND ELEMENTS OF TRUE ! ORBIT IF (STATE.EQ.'ABC') THEN CALL CALCTI ELSE CALL LOCATE (0,2) CALL WRITE_STRING('Specify all parameters first') END IF C GO TO 60 ! RETURN TO PGCURSE LINE C------------------------------------------------------------------ ELSE IF (COMND.EQ.'U') THEN ! UNDELETE A POINT DISTANCE = 1.E9 DO I = 1, NPTS DIST = (RA(I) - XR)**2 + (DEC(I) - YR)**2 IF (DIST.LT.DISTANCE) THEN DISTANCE = DIST II = I END IF END DO DELETE(II) = .FALSE. C------------------------------------------------------------------ ELSE IF (COMND.EQ.'V') THEN ! CHANGE RA, DEC MAX AND MIN ! IN PLOT CALL CHNGDISP C------------------------------------------------------------------ ELSE IF (COMND.EQ.'W') THEN ! WRITE EXPECTED SEPARATIONS, ! ETC., TO OUTFILE CALL WRITEMDL C------------------------------------------------------------------ ELSE IF (COMND.EQ.'X') THEN ! QUIT PROGRAM ALTOGETHER GO TO 80 C------------------------------------------------------------------ ELSE ! DISPLAY HELP COMND = '0' CALL LOCATE (0, 1) CALL WRITE_STRING ( '? help ' ) CALL LOCATE (0, 2) CALL WRITE_STRING ( 'A set major axis point' ) CALL LOCATE (0, 3) CALL WRITE_STRING ( 'B set minor axis point' ) CALL LOCATE (0, 4) CALL WRITE_STRING ( 'C set Center' ) CALL LOCATE (0, 5) CALL WRITE_STRING ( 'D Delete point' ) CALL LOCATE (0, 6) CALL WRITE_STRING ( 'E Error est. for true orbit' ) CALL LOCATE (0, 7) CALL WRITE_STRING ( 'F Fit apparent ellipse' ) CALL LOCATE (0, 8) CALL WRITE_STRING ( 'G local Grid search for chi^ min' ) CALL LOCATE (0, 9) CALL WRITE_STRING ( 'H Hardcopy' ) CALL LOCATE (0,10) CALL WRITE_STRING ( 'I Identify point' ) CALL LOCATE (0,11) CALL WRITE_STRING ( 'K Chi squared for given elements' ) CALL LOCATE (0,12) CALL WRITE_STRING ( 'L Least-squares true orbit' ) CALL LOCATE (0,13) CALL WRITE_STRING ( 'R Rotate ellipse' ) CALL LOCATE (0,14) CALL WRITE_STRING ( 'S (un)Show deleted points' ) CALL LOCATE (0,15) CALL WRITE_STRING ( 'T calc. True orbit elements' ) CALL LOCATE (0,16) CALL WRITE_STRING ( 'U Undelete point' ) CALL LOCATE (0,17) CALL WRITE_STRING ( 'V change View (max/min, symbols &c.)' ) CALL LOCATE (0,18) CALL WRITE_STRING ( 'W Write model RA, Dec, &c. to file' ) CALL LOCATE (0,20) CALL WRITE_STRING ( 'X Quit program' ) CALL LOCATE (0,21) CALL WRITE_STRING ( 'Q Quit to KEYIN level' ) C------------------------------------------------------------------ END IF C------------------------------------------------------------------ C END COMMAND INTERPRETATION C------------------------------------------------------------------ 80 CONTINUE C------------------------------------------------------------------ C CALCULATE THE PARAMETERS OF THE APPARENT ELLIPSE FROM THE C COORDINATES OF THE ENDS OF THE AXES AND OF THE CENTER IF (STATE.EQ.'AC') THEN D WRITE(OUTC,'(A,A)') ' After stmt 80: state = ',STATE A_APP = 0.5 * SQRT((ARA(2)-ARA(1))**2 + (ADEC(2)-ADEC(1))**2) A_ANGLE = ATAN2( (ARA(2)-ARA(1)) , (ADEC(2)-ADEC(1)) ) ELSE IF (STATE.EQ.'BC') THEN D WRITE(OUTC,'(A,A)') ' After stmt 80: state = ',STATE B_APP = 0.5 * SQRT((BRA(2)-BRA(1))**2 + (BDEC(2)-BDEC(1))**2) A_ANGLE = ATAN2( (BDEC(2)-BDEC(1)) , (BRA(2)-BRA(1)) ) ELSE IF (STATE.EQ.'ABC') THEN D WRITE(OUTC,'(A,A)') ' After stmt 80: state = ',STATE A_APP = 0.5 * SQRT((ARA(2)-ARA(1))**2 + (ADEC(2)-ADEC(1))**2) B_APP = 0.5 * SQRT((BRA(2)-BRA(1))**2 + (BDEC(2)-BDEC(1))**2) IF (B_APP.GT.A_APP) THEN ! IF B_APP > A_APP, SWAP AXES D WRITE(OUTC,*) ' B_APP > A_APP; swapping?!' DO I = 1,2 RADUMMY(I)= ARA(I) ARA(I) = BRA(I) BRA(I) = RADUMMY(I) DECDUMMY(I)= ADEC(I) ADEC(I) = BDEC(I) BDEC(I) = DECDUMMY(I) END DO RADUMMY(1)= A_APP A_APP = B_APP B_APP = RADUMMY(1) END IF A_ANGLE = ATAN2( (ARA(2)-ARA(1)) , (ADEC(2)-ADEC(1)) ) E_APP = SQRT( A_APP**2 - B_APP**2 ) / A_APP END IF IF (COMND.EQ.'Q'.OR.COMND.EQ.'X') THEN ! QUIT GO TO 100 ELSE IF (COMND.EQ.'0') THEN ! ASK BEFORE REFRESHING SCREEN CALL PGASK(.TRUE.) CALL PGADVANCE CALL PGASK(.FALSE.) GO TO 50 ELSE ! REFRESH SCREEN, LOOP BACK CALL PGASK(.FALSE.) CALL PGADVANCE GO TO 50 END IF 100 CONTINUE C CALL PGEND RETURN 900 WRITE(OUTC,501) ' Some kind of plot problem with PGBEGIN' RETURN 501 FORMAT(A) END SUBROUTINE CURPOS(OUTC,XR,YR) INTEGER*4 OUTC REAL*4 XR,YR WRITE(OUTC,*) ' Cursor position: ',XR,YR CALL PGASK(.TRUE.) CALL PGADVANCE CALL PGASK(.FALSE.) RETURN END