PROGRAM ELLIPSE C This program fits an ellipse to a set of points in polar coordinates. C It is intended for finding the apparent orbit of a binary. C Version 1.0: 19 Sept. 1989 Tom Armstrong C 1.1: 5 Jan. 1990 T.A. C Compile with: C del ellipse.obj C f77 ellipse.f -c -n4 -ep -w -s -lfgrex -lpg -lcgrex C del readpar.obj C f77 readpar.f -c -n4 -ep -w -s -lfgrex -lpg -lcgrex C del readdata.obj C f77 readdata.f -c -n4 -ep -w -s -lfgrex -lpg -lcgrex C del interpr.obj C f77 interpr.f -c -n4 -ep -w -s -lfgrex -lpg -lcgrex C del charnum.obj C f77 charnum.f -c -n4 -ep -w -s -lfgrex -lpg -lcgrex C del ctostr.obj C f77 ctostr.f -c -n4 -ep -w -s -lfgrex -lpg -lcgrex C del ctoi2.obj C f77 ctoi2.f -c -n4 -ep -w -s -lfgrex -lpg -lcgrex C del ctod2.obj C f77 ctod2.f -c -n4 -ep -w -s -lfgrex -lpg -lcgrex C del ctor2.obj C f77 ctor2.f -c -n4 -ep -w -s -lfgrex -lpg -lcgrex C del getpt.obj C f77 getpt.f -c -n4 -ep -w -s -lfgrex -lpg -lcgrex C del str$up.obj C f77 str$up.f -c -n4 -ep -w -s -lfgrex -lpg -lcgrex C del plotpars.obj C f77 plotpars.f -c -n4 -ep -w -s -lfgrex -lpg -lcgrex C del showdata.obj C f77 showdata.f -c -n4 -ep -w -s -lfgrex -lpg -lcgrex C del pointplt.obj C f77 pointplt.f -c -n4 -ep -w -s -lfgrex -lpg -lcgrex C del ellplot.obj C f77 ellplot.f -c -n4 -ep -w -s -lfgrex -lpg -lcgrex C del calcell.obj C f77 calcell.f -c -n4 -ep -w -s -lfgrex -lpg -lcgrex C del setell.obj C f77 setell.f -c -n4 -ep -w -s -lfgrex -lpg -lcgrex C del perpend.obj C f77 perpend.f -c -n4 -ep -w -s -lfgrex -lpg -lcgrex C del caps.obj C f77 caps.f -c -n4 -ep -w -s -lfgrex -lpg -lcgrex C del fitell.obj C f77 fitell.f -c -n4 -ep -w -s C del chisqr.obj C f77 chisqr.f -c -n4 -ep -w -s C del match.obj C f77 match.f -c -n4 -ep -w -s C del polint.obj C f77 polint.f -c -n4 -ep -w -s C del calcti.obj C f77 calcti.f -c -n4 -ep -w -s C del orbfit.obj C f77 orbfit.f -c -n4 -ep -w -s C del mrqmin.obj C f77 mrqmin.f -c -n4 -ep -w -s C del mrqcof.obj C f77 mrqcof.f -c -n4 -ep -w -s C del gaussj.obj C f77 gaussj.f -c -n4 -ep -w -s C del covsrt.obj C f77 covsrt.f -c -n4 -ep -w -s C del true2app.obj C f77 true2app.f -c -n4 -ep -w -s C del errest.obj C f77 errest.f -c -n4 -ep -w -s C Create library with: C 386LIB LIBELL.LIB -C ELLIPSE C 386LIB LIBELL.LIB -A READPAR C 386LIB LIBELL.LIB -A READDATA C 386LIB LIBELL.LIB -A INTERPR C 386LIB LIBELL.LIB -A CHARNUM C 386LIB LIBELL.LIB -A CTOSTR C 386LIB LIBELL.LIB -A CTOI2 C 386LIB LIBELL.LIB -A CTOD2 C 386LIB LIBELL.LIB -A CTOR2 C 386LIB LIBELL.LIB -A GETPT C 386LIB LIBELL.LIB -A STR$UP C 386LIB LIBELL.LIB -A PLOTPARS C 386LIB LIBELL.LIB -A SHOWDATA C 386LIB LIBELL.LIB -A POINTPLT C 386LIB LIBELL.LIB -A ELLPLOT C 386LIB LIBELL.LIB -A CALCELL C 386LIB LIBELL.LIB -A SETELL C 386LIB LIBELL.LIB -A PERPEND C 386LIB LIBELL.LIB -A CAPS C 386LIB LIBELL.LIB -A FITELL C 386LIB LIBELL.LIB -A MATCH C 386LIB LIBELL.LIB -A CHISQR C 386LIB LIBELL.LIB -A CALCTI C 386LIB LIBELL.LIB -A POLINT C 386LIB LIBELL.LIB -A INDEXX C 386LIB LIBELL.LIB -A ORBFIT C 386LIB LIBELL.LIB -A MRQMIN C 386LIB LIBELL.LIB -A MRQCOF C 386LIB LIBELL.LIB -A GAUSSJ C 386LIB LIBELL.LIB -A COVSRT C 386LIB LIBELL.LIB -A TRUE2APP C 386LIB LIBELL.LIB -A ERREST C COPY LIBELL.LIB L:\NDP\LIBELL.LIB C DELETE LIBELL.LIB C Link with: C f77 -o ellipse -n4 -lfgrex -lpg -lcgrex -lcit -lell tdate julda timec C error putout [all on one line] C PARAMETERS are in: PARAMETR.INC C DECLARATIONS are in: ELLIPSE.INC C All variables are in COMMON blocks LOGIC, CHAR, and NUM for logical, C character, and numeric variables. C SUBROUTINES in other files are: C READPAR.F CALLS KEYIN, READS IN PARAMETERS, OPENS FILES C READDATA.F READS IN DATA. CONTAINS SUBROUTINES READELL, C STR$UPCASE, GETPT, CTOR2, CTOD2, CTOI2, C CTOSTR,CHARNUM,INTERPR. ALSO CALLS C JULDA C SHOWDATA.F CALLED BY SETELL & SHOWFIT. CONTAINS: C PLOTPARS CALCULATES PLOT MINIMUM, MAXIMUM; ALSO C CALCULATES ERROR BAR COORDINATES C SHOWDATA DISPLAYS DATA AND UNCERTAINTIES C POINTPLT DOES ACTUAL POINT PLOTTING C PLOTELL PLOTS THE ELLIPSE C CALCELL CALCULATES THE ELLIPSE TO BE PLOTTED C SETELL.F SETS INITIAL ELLIPSE PARAMETERS INTERACTIVELY C FITELL.F FITS ELLIPSE BY GRADIENT SEARCH, GRID SEARCH. C ALSO INCLUDES: C INDEXX INDEXING PROG. FROM NUMERICAL RECIPES C SHOWFIT.F DISPLAYS DATA AND FITTED ELLIPSE C CALCTI.F CALCULATES THIELE-INNES CONSTS, ETC.; C WRITES FILES C ORBFIT.F DOES LEAST-SQUARES FIT TO TRUE ORBIT C MRQMIN.F CALLED BY ORBFIT TO SET UP LEAST-SQUARES C MRQCOF.F CALLED BY MRQMIN TO DO AN ITERATION OF LEAST- C SQUARES FIT C TRUE2APP.F CALCULATES APPARENT POSITIONS FOR GIVEN OBSER- C VATION TIMES FROM ELEMENTS OF TRUE ORBIT C ERREST.F ESTIMATES ERRORS IN LEAST-SQUARES FIT C LIST OF VARIABLES: C ELEMENTS OF ELLIPSE OF APPARENT ORBIT: C A_APP SEMI-MAJOR AXIS OF APPARENT ORBIT(mas) C ARA(2) ENDS OF MAJOR AXIS OF APPARENT ORBIT (from cursor) C ADEC(2) C B_APP SEMI-MINOR AXIS OF APPARENT ORBIT (mas) C BRA(2) ENDS OF MINOR AXIS OF APPARENT ORBIT (from cursor) C BDEC(2) C E_APP ELLIPTICITY C R_ELL(N) RADIAL COORDINATE IN ELLIPSE (mas) C THETA_ELL(N) ANGULAR COORDINATE IN ELLIPSE (deg) C RA_ELL(N) R.A. ALONG CALCULATED ELLIPSE C DEC_ELL(N) DEC. ALONG CALCULATED ELLIPSE C THETA_PT(N) POS. ANGLE ALONG CALCULATED ELLIPSE (rad) C THETA_LINE POS. ANGLE OF CTR OF CALC. ELLIPSE (rad) C A_ANGLE POSITION ANGLE OF MAJ. AX. (deg in KEYIN, else rad) C CRA X-COORDINATE OF ELLIPSE CENTER (mas) C CDEC Y-COORDINATE OF ELLIPSE CENTER (mas) C RAINT USED IN MINOR AXIS CALCULATION C DECINT USED IN MINOR AXIS CALCULATION C DIST1, DIST2 USED IN MAJOR AXIS CALCULATION C AREA_TIME AREA ON SKY SWEPT OUT BY LINE CONNECTING C STARS (mas^2/day) C ELEMENTS OF TRUE ORBIT: C A_TRUE SEMI-MAJOR AXIS (mas) C E_TRUE ECCENTRICITY C PERIOD IN DAYS C EPOCH TIME OF PERIASTRON PASSAGE (days) C INCLIN INCLINATION ANGLE (deg in KEYIN, else rad) C ASC_NODE POS. ANGLE OF ASCENDING NODE (deg in KEYIN, else rad) C ARG_PERI ARGUMENT OF PERIASTRON (deg in KEYIN, else rad) C A(7) ARRAY OF SEVEN ELEMENTS, USED BY MRQMIN ETC. C M_ANOM MEAN ANOMALY (radians) C E_ANOM ECCENTRIC ANOMALY (radians) C T_ANOM TRUE ANOMALY (radians) C R_TRUE SEPARATION IN TRUE ORBIT (mas) C DATA: C ROBS(N) OBSERVED RADIAL COORDINATE (mas) C DROBS(N) UNCERTAINTY IN RADIAL COORDINATE OF OBSERVATION C THETAOBS(N) OBSERVED ANGULAR COORDINATE (deg; converted C to radians in ) C DTHETAOBS(N) UNCERTAINTY IN ANGULAR COORD. OF OBSERVATION C (deg; converted to radians in ) C TSTRING(N) TIME IN yymmdd FORMAT C UT(N) UT OF OBSERVATION (IN SECONDS!) C JULDAY(N) JULIAN DAY OF OBSERVATION (from JULDA) C GMST(N) GMST OF OBSERVATION (from JULDA) C TOBS(N) TIME OF OBSERVATION (days from Jan. 1 of C first year of observations) C TOBSREF TIME OF REFERENCE OBSERVATION (for time calc.) C DELETE(N) IS Nth POINT DELETED? C FIRSTYR YEAR OF FIRST OBSERVATION C NPTS NUMBER OF POINTS IN DATA SET C NFIT NUMBER OF POINTS TO FIT ELLIPSE TO C STARNAME C FK5 STAR'S FK5 NUMBER C T_ZERO TIME OF PERIASTRON (years from Jan. 1) C ORBSENSE "Retrograde" (caps or lc) => RETROGRADE ORBIT; C OTHERWISE PROGRADE C SENSE -1 => RETROGRADE; +1 => PROGRADE C RA_EST(N) PREDICTED RA, DEC FOR TIMES AT WHICH THE DATA C DEC_EST(N) WERE TAKEN C R_EST PREDICTED RADIAL COORDINATES FOR TIMES AT C THETA_EST WHICH THE DATA WERE TAKEN C THIELE-INNES CONSTANTS C X_TH_1 X-COORDINATE OF PERIASTRON C Y_TH_1 Y-COORDINATE OF PERIASTRON C X_TH_2 X-COORDINATE OF V=90 deg POSITION C Y_TH_2 Y-COORDINATE OF V=90 deg POSITION C DISPLAY: C RA(N), DEC(N) COORDINATES OF OBSERVED POINTS (mas) C RADR, DECDR COORDS OF ENDPOINTS OF RADIAL UNCERT C RADTH, DECDTH COORDS OF ENDPOINTS OF ANGULAR UNCERT C RAMIN, RAMAX MINIMUM, MAXIMUM ALONG RA AXIS C DECMIN, DECMAX MINIMUM, MAXIMUM ALONG DEC AXIS C RALO, RAHI MIN, MAX OF GRAPH C DECLO, DECHI MIN, MAX OF GRAPH C NFONT FONT C LINEWD C PLOTFILE PLOT FILE C XR, YR CURSOR POSITION C TITLE(5) C COMND COMMAND CHARACTER RETURNED BY PGCURSE C STATE = 'N' IF NO POINTS HAVE BEEN SET C = 'A', 'B', 'C' IF ONLY ONE POINT HAS BEEN SET C = 'AB', 'AC', 'BC' IF TWO POINTS HAVE BEEN SET C = 'ABC' IF THREE POINTS HAVE BEEN SET SUCH THAT C BOTH AXES AND THE CENTER ARE DETERMINED C SHOWDEL .TRUE. IF DELETED POINTS SHOWN C ELLFIT .TRUE. IF FIT HAS BEEN DONE; DISPLAY CALCULATED C POSITIONS TO COMPARE WITH OBSERVED POSITIONS C FILES: C DATAFILE INPUT DATA FILE C OUTFILE OUTPUT FILE WITH T-H CONSTANTS, ETC. C MISCELLANEOUS: C INCLUDE 'PARAMETR.INC' INCLUDE 'ELLIPSE.INC' C DATA PARS/ 'RAmax', 'DECmax', 'MArk', 'DOts', 'LInewidth', 1 'Font', 'PLot_dev', 11*' ', 'INPutfile', 11*' ', 2 'OUTputfile', 11*' ', 3 'AAParent','BAPparent','EAPparent', 4 'AANgle','CRA','CDEC', 5 'EPoch','INClination','ARGperi','ASCending', 6 'ATrue','ETrue','PERiod' 7 / DATA VALS/ 0D0, 0D0, -1D0, 0D0, 1D0, 1 1D0, '/LPM', 11*' ','?',11*' ', 2 '?',11*' ', 3 13*-1000D0 4 / DATA PLOTFILE /'PLOT.NNN'/ DATA ENDMRK/ '/' / DATA BLANK/' '/ DATA MINUS/'-'/ C----------------------------------------------------------------------- C Read file names and parameters: C----------------------------------------------------------------------- 10 CALL READPAR C---------------------------------------------------------------------- C Read data file C Data needed are: STARNAME and PERIOD (first record); then C separation, position angle, time for each point, and C uncertainties. C Calculate the RA, DEC coordinates of the points and of the C endpoints of the error bars: C---------------------------------------------------------------------- CALL READDATA C---------------------------------------------------------------------- C Calculate minimum and maximum for display, and the coordinates C of the error bars: C---------------------------------------------------------------------- CALL PLOTPARS C---------------------------------------------------------------------- C Display the data and, if available, the apparent and true orbits. C SETELL acts as a shell for interpreting various commands. They C include: C Get the initial parameters of the ellipse interactively. Model C parameters needed are: C major axis, eccentricity, position angle of major axis, C coordinates of center C Calculate the true orbit by the Thiele-Innes method C The least-squares true orbit fit is also run from here, using either C the Thiele-Innes constants or values from KEYIN. C---------------------------------------------------------------------- CALL SETELL C---------------------------------------------------------------------- C On return, if COMND = Q, set the KEYIN parameters from the inter- C actively set apparent ellipse and return to READPAR. If COMND = X, C exit from the program entirely. C---------------------------------------------------------------------- IF (COMND.EQ.'Q') THEN IF (STATE.EQ.'AC'.OR.STATE.EQ.'ABC') THEN VALS(43) = A_APP END IF IF (STATE.EQ.'BC'.OR.STATE.EQ.'ABC') THEN VALS(44) = B_APP END IF IF (STATE.EQ.'ABC') THEN VALS(45) = E_APP END IF IF (STATE.EQ.'AC'.OR.STATE.EQ.'BC'.OR.STATE.EQ.'ABC') THEN VALS(46) = A_ANGLE / RPDEG VALS(47) = CRA VALS(48) = CDEC END IF GO TO 10 END IF C---------------------------------------------------------------------- GO TO 3000 C---------------------------------------------------------------------- C Error exits C---------------------------------------------------------------------- 900 CALL ERROR('Cannot open PLOTFILE') 3000 CALL PGEND END