C	NEOWISER-when.f - list when a given position will be observed by WISE
C	g77 -o NEOWISER-when NEOWISER-when.f ~/Ned.a
C	updated 14-MAR-2015 for 2015 scan parameters
C
	IMPLICIT NONE
	REAL*4 E(3),C(3),A(3,3),SC(3)
	REAL*8 T01JAN14,T
	LOGICAL CL
	CHARACTER*255 INPUT,ARG
	INTEGER MON,IH,IDAY,IY,IM
	REAL S,ELONG,DTR
	REAL BIAS
	REAL DIHEDRAL
	INTEGER NTOKEN,IARGC
	INTEGER I,K,IDL
	REAL TH,SL
	INTEGER NC,MC
	REAL SUNLONG
C
	DTR = 45/ATAN(1.)
	CALL TIME_TAG(2014,1,1,0,0,0.,T01JAN14)
	NTOKEN = IARGC()
	CL = NTOKEN.GT.0
C
	IF (CL) THEN
          DO I=1,NTOKEN
            CALL GETARG(I,ARG)
            MC = LNBLNK(ARG)
	    CALL CONVERT_POSITION(ARG(1:MC),C)
	    CALL CTOE(C,E)
	    ELONG = DTR*ATAN2(E(2),E(1))
C
	    DO K=-15,1096
	      T = T01JAN14+K*864.D2 + 100
	      CALL T68_TO_UTC(T,IY,MON,IDAY,IH,IM,S)
	      BIAS = 0
	      DIHEDRAL = 2.5
	      IF (K.GT.385) THEN
	        BIAS = 4.13
	        DIHEDRAL = 4.13
	      ENDIF
	      TH = 360*K/365.2425
	      SL = -79.596+1.914*SIN(TH/DTR)-0.088*COS(TH/DTR)+TH
	      DO IDL=-1,1,2
C 2014:	        SUNLONG = ELONG+92.5*IDL
C 2015: ELONG = SUNLONG+98.26 or SUNLONG-90
	        SUNLONG = ELONG+(90+DIHEDRAL)*IDL-BIAS
	        IF (ABS(MOD(SL-SUNLONG+540,360.0)-180).LT.0.5) THEN
	          WRITE (*,901) MON,IDAY,IY
	        ENDIF
	      ENDDO
	    ENDDO
C
	  ENDDO
	  STOP
        ENDIF
10	WRITE (*,900) 'Enter Type(E/C/G),long,lat:'
900     FORMAT(A,$)
	  READ (*,FMT='(A)',END=100) INPUT
	  NC = LNBLNK(INPUT)
	  CALL CONVERT_POSITION(INPUT(1:NC),C)
	  CALL CTOE(C,E)
	  ELONG = DTR*ATAN2(E(2),E(1))
C
	  DO K=-15,1096
	    T = T01JAN14+K*864.D2 + 100
	    CALL T68_TO_UTC(T,IY,MON,IDAY,IH,IM,S)
	    BIAS = 0
	    DIHEDRAL = 2.5
	    IF (K.GT.385) THEN
	      BIAS = 4.13
	      DIHEDRAL = 4.13
	    ENDIF
	    TH = 360*K/365.2425
	    SL = -79.596+1.914*SIN(TH/DTR)-0.088*COS(TH/DTR)+TH
	    DO IDL=-1,1,2
C 2014:	        SUNLONG = ELONG+92.5*IDL
C 2015: ELONG = SUNLONG+98.26 or SUNLONG-90
	      SUNLONG = ELONG+(90+DIHEDRAL)*IDL-BIAS
	      IF (ABS(MOD(SL-SUNLONG+540,360.0)-180).LT.0.5) THEN
	        WRITE (*,901) MON,IDAY,IY
	      ENDIF
	    ENDDO
	  ENDDO
C
	GO TO 10
901	FORMAT(I2,1H/,I2,1H/,I4,I3,1H:,I2,1H:,I2)
100	WRITE (*,FMT='(A)')
	STOP
	END

	SUBROUTINE CONVERT_POSITION(INPUT,C)
C	converts a string into a position
	REAL*4 G(3),E(3),C(3),LNG,LAT,PMT(3,3)
	CHARACTER*1 SW
	CHARACTER*(*) INPUT
	CHARACTER*17 STR
	DTR = 45./ATAN(1.)
	PI = 180/DTR
901	FORMAT('(l,b) = ',2F8.3)
902	FORMAT(A)
903	FORMAT('(lambda,beta) = ',2F8.3)
904	FORMAT(A,3F10.1)
	READ (INPUT,904) SW,LNG,LAT,EPOCH
	IF (SW.EQ.'E') THEN
	  E(3) = SIN(LAT/DTR)
	  E(1) = COS(LAT/DTR)*COS(LNG/DTR)
	  E(2) = COS(LAT/DTR)*SIN(LNG/DTR)
	  WRITE (*,903) LNG,LAT
	  CALL ETOC(E,C)
	  CALL CTOG(C,G)
	  CALL C_TO_STRING(C,0,STR)
	  LNG = DTR*ATAN2(G(2),G(1))
	  IF (LNG.LT.0.) LNG = LNG+360
	  LAT = DTR*ATAN2(G(3),SQRT(G(1)**2+G(2)**2))
	  WRITE (*,901) LNG,LAT
	  WRITE (*,902) STR
	ELSE IF (SW.EQ.'C' .OR. SW.EQ.'D') THEN
	  IF (SW.EQ.'C') THEN
	    LAT = DMS(LAT)
	    LNG = 15.*DMS(LNG)
	  ELSE
	    LAT = LAT/DTR
	    LNG = LNG/DTR
	  ENDIF
	  C(3) = SIN(LAT)
	  C(1) = COS(LAT)*COS(LNG)
	  C(2) = COS(LAT)*SIN(LNG)
	  IF (EPOCH.NE.0.0 .AND. ABS(EPOCH-2000.).GT.0.1) THEN
	    CALL PREMAT(EPOCH,PMT)
	    CALL VMAT(PMT,C,E)
	    DO I=1,3
	      C(I) = C(I)+E(I)
	    ENDDO
	  ENDIF
	  CALL CTOE(C,E)
	  CALL CTOG(C,G)
	  LNG = DTR*ATAN2(G(2),G(1))
	  IF (LNG.LT.0.) LNG = LNG+360
	  LAT = DTR*ATAN2(G(3),SQRT(G(1)**2+G(2)**2))
	  WRITE (*,901) LNG,LAT
	  CALL C_TO_STRING(C,0,STR)
	  WRITE (*,902) STR
	  LNG = DTR*ATAN2(E(2),E(1))
	  IF (LNG.LT.0.) LNG = LNG+360
	  LAT = DTR*ATAN2(E(3),SQRT(E(1)**2+E(2)**2))
	  WRITE (*,903) LNG,LAT
	ELSE IF (SW.EQ.'G') THEN
	  G(3) = SIN(LAT/DTR)
	  G(1) = COS(LAT/DTR)*COS(LNG/DTR)
	  G(2) = COS(LAT/DTR)*SIN(LNG/DTR)
	  CALL GTOC(G,C)
	  CALL CTOE(C,E)
	  WRITE (*,901) LNG,LAT
	  CALL C_TO_STRING(C,0,STR)
	  WRITE (*,902) STR
	  LNG = DTR*ATAN2(E(2),E(1))
	  IF (LNG.LT.0.) LNG = LNG+360
	  LAT = DTR*ATAN2(E(3),SQRT(E(1)**2+E(2)**2))
	  WRITE (*,903) LNG,LAT
	ENDIF
	RETURN
100	STOP
	END

C	routines from nedlib.f below here:

	SUBROUTINE CTOE(C,E)
	REAL*4 C(3),E(3),TEMP
C
C	C	INPUT	vector in celestial coordinates
C	E	OUTPUT	vector in ecliptic coordinates
C
	PARAMETER (COBL=.917482058,SOBL=.39777716)
	E(1)=C(1)
	TEMP=COBL*C(3)-SOBL*C(2)
	E(2)=COBL*C(2)+SOBL*C(3)
	E(3)=TEMP
	RETURN
	END

	SUBROUTINE ETOC(E,C)
	REAL*4 E(3),C(3),TEMP
C
C	E	INPUT	vector in ecliptic coordinates (J2000.0)
C	C	OUTPUT	vector in celestial coordinates (J2000.0)
C
	PARAMETER (COBL=.917482058,SOBL=.39777716)
	C(1)=E(1)
	TEMP=COBL*E(3)+SOBL*E(2)
	C(2)=COBL*E(2)-SOBL*E(3)
	C(3)=TEMP
	RETURN
	END

	SUBROUTINE GTOC(G,C2000)
C
C	CONVERTS VECTOR G FROM GALACTIC COORDINATES to J2000 in C
C
	DIMENSION G(3),C2000(3),R(3,3)
	DATA R/ -0.0548809,  0.4941080, -0.8676667,
	1       -0.8734369, -0.4448323, -0.1980717,
 	2       -0.4838349,  0.7469817,  0.4559848/
C
	CALL VMAT(R,G,C2000)
	RETURN
	END

	SUBROUTINE CTOG(C2000,G)
C
C	CONVERTS VECTOR C FROM CELESTIAL COORDINATES (2000.0)
C	TO GALACTIC COORDINATE VECTOR G
C
	DIMENSION G(3),C2000(3),R(3,3)
	DATA R/ -0.0548809,  0.4941080, -0.8676667,
	1       -0.8734369, -0.4448323, -0.1980717,
 	2       -0.4838349,  0.7469817,  0.4559848/
	CALL MAV(R,C2000,G)
	RETURN
	END

	REAL*4 FUNCTION DELTA_T(T68)
C
C	T68 is REAL*8 seconds since 24 May 1968, UT 00:00:00
C	which was JD 2440000.5
C
C	DELTA_T is the difference between IAT (or ET or TDT) and UTC
C	since 5/24/68
C
C	the value returned includes the 32.18 seconds between terrestrial
C	time and TAI
C
C	The seconds to be used are ephemeris seconds, or IAT seconds
C	Therefore the conversion from UT date and time must allow for the
C	rates adjustments used before 1972 and the leap seconds used later
C	This routine first calculates UTC ignoring these corrections, then
C	looks in a table to interpolate for DT = ET - UT
C	Linear interpolation is used, so two entries are made for each
C	leap second.  Entry for JDATE should be the JD at noon the day 
C	after the leap second.
C
	IMPLICIT NONE
	INTEGER*4 NPT
	PARAMETER (NPT=32)
	INTEGER*4 JDATE(NPT),I,NL,NH,NTRY
	REAL*8 DT(NPT),DTNOW,SLOPE
	LOGICAL SAVED/.FALSE./
	INTEGER*4 JD68/2440001/
	REAL*8 DT68/38.66/
	REAL*8 T68,TAB68(NPT)
C
	SAVE SAVED,TAB68
C
	DATA (JDATE(I),DT(I),I=1,NPT)/
	1		2440001,	38.66,	      	!	5/24/1968
	1		2445517,	53.18,		!	7/1/1983
	1		2445517,	54.18,		!	7/1/1983
	1		2446248,	54.18,		!	7/1/1985
	1		2446248,	55.18,		!	7/1/1985
	1		2447162,	55.18,		!       1/1/1988
	1		2447162,	56.18,		!	1/1/1988
	1		2447893,	56.18,		!	1/1/1990
	1		2447893,	57.18,		!	1/1/1990
	1		2448258,	57.18,		!	1/1/1991
	1		2448258,	58.18,		!	1/1/1991
	1		2448805,	58.18,		!	7/1/1992
	1		2448805,	59.18,		!	7/1/1992
	1		2449170,	59.18,		!     	7/1/1993
	1		2449170,	60.18,		!     	7/1/1993
	1		2449535,	60.18,		!	7/1/1994
	1		2449535,	61.18,		!	7/1/1994
	1		2450084,	61.18,		!	1/1/1996
	1		2450084,	62.18,		!	1/1/1996
	1		2450631,	62.18,		!	7/1/1997
	1		2450631,	63.18,		!	7/1/1997
	1		2451180,	63.18,		!	1/1/1999
	1		2451180,	64.18,		!	1/1/1999
	1		2453737,	64.18,		!	1/1/2006
	1		2453737,	65.18,		!	1/1/2006
	1		2454833,	65.18,		!	1/1/2009
	1		2454833,	66.18,		!	1/1/2009
	1		2456110,	66.18,		!	7/1/2012
	1		2456110,	67.18,		!	7/1/2012
	1		2457205,	67.18,		!	7/1/2015
	1		2457205,	68.18,		!	7/1/2015
	1		2488070,	68.18/		!	1/1/2100
C
C---------------------------------------------------
	IF (.NOT. SAVED) THEN
C				 YYDDDHHMMSSFFF
	  SAVED = .TRUE.
	  DO I=1,NPT
	    TAB68(I)=(JDATE(I)-JD68)*864.D2 + (DT(I)-DT68)
	  ENDDO
	ENDIF
C
C	Find position in Table
C
	NL=1
	NH=NPT-1
	DO WHILE (NL.LT.NH)
	  NTRY=(NL+NH+1)/2
	  IF(TAB68(NTRY).GT.T68) THEN
	    NH=NTRY-1
	  ELSE
	    NL=NTRY
	  ENDIF
	ENDDO
	SLOPE = (DT(NL+1)-DT(NL))/(TAB68(NL+1)-TAB68(NL))
	DTNOW=DT(NL)+ (T68-TAB68(NL)) * SLOPE
	DELTA_T = (DTNOW-DT68)	! time correction ET-UTC in seconds
	RETURN
	END

	SUBROUTINE TIME_TAG(IY,IM,ID,IH,MINUTE,SECOND,T)
C
C	CONVERTS YEAR, MONTH, DAY, HOUR, MINUTE, SECOND INTO
C	REAL*8 TIMETAG T, SECONDS SINCE 5/24/68 00:00:00 UT
C	The time tag T is REAL*8 seconds since 24 May 1968, UT 00:00:00
C	which was JD 2440000.5
C
C	The seconds to be used are ephemris seconds, or IAT seconds
C	Therefore the conversion from UT date and time must allow for the
C	rates adjustments used before 19?? and the leap seconds used later
C	This routine first calculates UTC ignoring these corrections, then
C	looks in a table to interpolate for DT = ET - UT
C	Linear interpolation is used, so two entries are made for each
C	leap second.  Entry for JDATE should be the JD at noon the day 
C	after the leap second.
C
	REAL*8 T,DT
	REAL*4 SECOND,DELTA_T
	INTEGER*4 JD,I367,JULIAN,JDZERO
	DATA JDZERO/2440001/			!JD at 5/24/68 UT noon
	DATA I367/367/
	JD=-7*(IY+(IM+9)/12)/4-3*((IY+(IM-9)/7)/100+1)/4+275*IM/9+ID
	JULIAN=JD+I367*IY+1721029
C	lookup DT at UT noon to avoid leap seconds
	T=(JULIAN-JDZERO)*864.D2+432.D2
	DT = DELTA_T(T)
	T=(JULIAN-JDZERO)*864.D2+IH*36.D2+MINUTE*6.D1+SECOND+DT
	RETURN
	END

	SUBROUTINE T68_TO_UTC(T,IY,MON,ID,IH,MINUTE,SECOND)
C
C	CONVERTS TIME TAG T INTO UT TIME IH:MINUTE:SECOND ON IM/ID/IY
C	T is seconds since 5/24/68 allowing for leap seconds
C
	REAL*8 T,TT,DT,DTT
	INTEGER*4 MAXDAY(12)
	REAL*8 DAY
	DATA MAXDAY/31,28,31,30,31,30,31,31,30,31,30,31/
	DAY = T/864.D2+140
	IY = DAY/365.25D0+1968
	DAY = MOD(DAY,365.25D0)
	MON = DAY/30.5D0+1.
	ID = MOD(DAY,30.5D0)+1.
10	MAXD = MAXDAY(MON)
C	check for Leap year Februaries:
	IF (MON.EQ.2 .AND. ((MOD(IY,4).EQ.0 .AND. MOD(IY,100).NE.0) .OR.
	1	MOD(IY,400).EQ.0) ) MAXD=29
	IF (ID.GT.MAXD) THEN
	  ID = ID-MAXD
	  MON = MON+1
	  IF (MON.GT.12) THEN
	    MON = 1
	    IY = IY+1
	  ENDIF
	  GO TO 10
	ENDIF
	CALL TIME_TAG(IY,MON,ID,0,0,1.,TT)
	DT = T-TT+1.D0
C	MAKE SURE THE DAY IS RIGHT
	IF (DT.LT.86400.D0) THEN
	  IH = DT/36.D2
	  DT = DT-36.D2*IH
	  MINUTE = DT/60.D0
	  SECOND = DT-60.D0*MINUTE
	  RETURN
	ENDIF
C	CHECK FOR A LEAP SECOND
	IF (DT.LT.86401.D0) THEN
	  CALL TIME_TAG(IY,MON,ID+1,0,0,1.,TT)
	  DTT = T-TT+1.D0
	  IF (DTT.GE.0.) THEN
	    ID = ID+1
	    GO TO 10
	  ENDIF
	  IH = 23
	  MINUTE = 59
	  SECOND = DT - 86340.D0
	  RETURN
	ENDIF
	ID = ID+DT/86401.D0
	GO TO 10
	END

	REAL*4 FUNCTION DMS(R)
C
C	CONVERTS INPUT REAL*4 R FROM "DD.MMSS" INTO RADIANS
C
	REAL*4 R
	PARAMETER (PI=3.1415927,TWOPI=6.2831853,DTR=57.29577951)
	V=ABS(R)
	ID=V
	V=100.*(V-FLOAT(ID))
	IM=V
	V=100.*(V-FLOAT(IM))
	V=FLOAT(ID)+FLOAT(IM)/60.+V/3600.
	IF(R.LT.0.) V=-V
	DMS=V/DTR
	RETURN
	END

	SUBROUTINE C_TO_STRING(DCEL,NDIG,STRING)
C
C	CONVERTS RA+DEC UNIT VECTOR DCEL INTO A FORMATTED STRING FOR OUTPUT
C
C	DCEL IS REAL*4 (3) INPUT ONLY
C	NDIG IS AN INTEGER INPUT THAT CONTROLS PRECISION
C	FOR NDIG = 0, THE OUTPUT IS HHhMMmSS -DDdMM.M = 17 characters
C	           1,               HHhMMmSS.S -DDdMM'SS = 20 chars
C		   2,		    HHhMMmSS.SS -DDdMM'SS.S = 23 chars
C
C	ROUNDING CAN GIVE 60.0 for the seconds PART
C
	CHARACTER*(*) STRING
	CHARACTER*1 IDS
	PARAMETER (PI=3.1415927,TWOPI=6.2831853,DTR=57.295778)
	REAL*4 DCEL(3)
	RA=ATAN2(DCEL(2),DCEL(1))
	IF(RA.LT.0.) RA=RA+TWOPI
	DEC=ATAN2(DCEL(3),SQRT(DCEL(1)**2+DCEL(2)**2))
	IF(DEC.LT.0.) THEN
	  IDS = '-'
	  DEC=-DEC
	ELSE
	  IDS = ' '
	ENDIF
	CALL R TO DMS(RA/15.,IH,IRM,SR)
	CALL R TO DMS(DEC,ID,IDM,SD)
	IF (NDIG.EQ.2) THEN
	  WRITE (STRING,900) IH,IRM,SR,IDS,ID,IDM,SD
	  IS=13
900	FORMAT(I2,1Hh,I2,1Hm,F5.2,1X,A,I2,1Hd,I2,1H',F4.1)
	ELSE IF (NDIG.EQ.1) THEN
	  ISD=SD+.5
	  IF(ISD.EQ.60) THEN
	    ISD=0
	    IDM=IDM+1
	    IF(IDM.EQ.60) THEN
	      IDM=0
	      ID=ID+1
	    ENDIF
	  ENDIF
	  WRITE (STRING,901) IH,IRM,SR,IDS,ID,IDM,ISD
	  IS=12
901	FORMAT(I2,1Hh,I2,1Hm,F4.1,1X,A,I2,1Hd,I2,1H',I2)
	ELSE
	  RM=IDM+SD/60.
	  ISR=SR+.5
	  IF(ISR.EQ.60) THEN
	    ISR=0
	    IRM=IRM+1
	    IF(IRM.EQ.60) THEN
	      IRM=0
	      IH=IH+1
	      IF(IH.EQ.24) IH=0
	    ENDIF
	  ENDIF
	  WRITE (STRING,902) IH,IRM,ISR,IDS,ID,RM
	  IS=10
902	FORMAT(I2,1Hh,I2,1Hm,I2,1X,A,I2,1Hd,F4.1)
	ENDIF
	IF(STRING(IS:IS+1).EQ.'- ') STRING(IS:IS+1) = ' -'
	RETURN
	END

	SUBROUTINE PREMAT(EPOCH,R)
C
C	COMPUTES MATRIX FOR PRECESSION FROM EPOCH TO 2000.
C
C	V(EPOCH)=(I+R)V(2000)
C
C	             T
C	V(2000)=(I+R) V(EPOCH)
C
	REAL*4 R(3,3),RE(3,3),R2(3,3)
	INTEGER IFF
	DATA IFF /0/
	SAVE IFF,R2
	IF (IFF.EQ.0) THEN
	  IFF = 1
	  CALL PREMAT50(2000.0,R2)
	ENDIF
	CALL PREMAT50(EPOCH,RE)
C						      T
C	V(EPOCH)=(I+R)V(2000)=(I+RE)V(50)=(I+RE)(I+R2) V(2000)
C		      T	       T
C	SO R = RE + R2  + RE*R2
C
	DO 30 J=1,3
	  DO 20 I=1,3
	    R(I,J)=RE(I,J)+R2(J,I)
	    DO 10 K=1,3
	      R(I,J) = R(I,J) + RE(I,K)*R2(J,K)
10	    CONTINUE
20	  CONTINUE
30	CONTINUE
	RETURN
	END

	SUBROUTINE PREMAT50(EPOCH,R)
C
C	COMPUTES MATRIX FOR PRECESSION FROM EPOCH TO 1950.
C
C	V(EPOCH)=(I+R)V(50)
C
C	           T
C	V(50)=(I+R) V(EPOCH)
C
	REAL*4 R(3,3),AP(3,3),BP(3,3),CP(3,3),P(3,3)
	DTR=45./ATAN(1.)
	T=(EPOCH-1950.0)/100.
	ZETA=T*(.6402633+T*(.0000839+T*.000005))/DTR
	Z=ZETA+.0002197*T**2/DTR
	THETA=T*(0.5567376-T*(.0001183+T*.0000117))/DTR
	DO 10 I=1,3
	AP(I,3)=0.
	AP(3,I)=0.
	CP(I,3)=0. 
	CP(3,I)=0.
	BP(I,2)=0.
10	BP(2,I)=0.
	AP(1,1)=-2.*SIN(Z/2.)**2
C
C	COS(Z)-1 = -2*SIN(Z/2)**2
C
	AP(2,2)=AP(1,1)
	AP(1,2)=-SIN(Z)
	AP(2,1)=-AP(1,2)
	BP(1,1)=-2.*SIN(THETA/2.)**2
	BP(3,3)=BP(1,1)
	BP(1,3)=-SIN(THETA)
	BP(3,1)=-BP(1,3)
	CP(1,1)=-2.*SIN(ZETA/2.)**2
	CP(2,2)=CP(1,1)
	CP(1,2)=-SIN(ZETA)
	CP(2,1)=-CP(1,2)
C
C	AP=A-I, BP=B-I, CP=C-I
C	NOW COMPUTE BC-I=AP*BP+AP+BP
C
	DO 30 I=1,3
	DO 30 J=1,3
	S=0.
	DO 20 K=1,3
20	S=S+BP(I,K)*CP(K,J)
30	P(I,J)=S+BP(I,J)+CP(I,J)
C
C	NOW COMPUTE R=ABC-I=AP*(BC-I)+AP+(BC-I)
C
	DO 50 I=1,3
	DO 50 J=1,3
	S=0.
	DO 40 K=1,3
40	S=S+AP(I,K)*P(K,J)
50	R(I,J)=S+AP(I,J)+P(I,J)
	RETURN
	END


	SUBROUTINE MAV(A,VSC,VCEL)
C
C	MULTIPLIES MATRIX A TIMES VECTOR VSC GIVING VECTOR VCEL
C
	REAL*4 A(3,3),VSC(3),VCEL(3),SC(3),VC
	DO 5 I=1,3
5	SC(I)=VSC(I)
	DO 15 I=1,3
	VC=0.
	DO 10 J=1,3
10	VC=VC+A(I,J)*SC(J)
15	VCEL(I)=VC
	RETURN
	END

	SUBROUTINE VMAT(A,VCEL,VSC)
C
C	MULTIPLIES TRANSPOSE OF A TIMES VCEL GIVING VSC
C
	REAL*4 A(3,3),VCEL(3),VSC(3),CEL(3),SC
	DO 5 I=1,3
5	CEL(I)=VCEL(I)
	DO 15 I=1,3
	SC=0.
	DO 10 J=1,3
10	SC=SC+A(J,I)*CEL(J)
15	VSC(I)=SC
	RETURN
	END

	SUBROUTINE RTODMS(R,ID,IM,S)
C
C	CONVERTS INPUT RADIANS R INTO INTEGER DEGREES, MINUTES AND REAL 
C	SECONDS IN ID,IM,S.  FAILS FOR -0 30' ETC BECAUSE VAX DOESN'T 
C	HAVE A MINUS ZERO
C
	PARAMETER (PI=3.1415927,TWOPI=6.2831853,DTR=57.29577951)
	V=ABS(R)*DTR
	ID=V
	V=60.*(V-FLOAT(ID))
	IM=V
	S=60.*(V-FLOAT(IM))
	IF(R.LT.0.) ID=-ID
	RETURN
	END