SUBROUTINE SRPUP(BETAPRIME, ALPHA, ATIME, MASS, AREA, CR, 
     *		ORB, DT)

	DOUBLE PRECISION BETAPRIME, ALPHA, MASS, ORB(6), DT, AREA, 
     *		CR, DEDT, DWDT, S, F, PHASE, CMIN, DELTA,
     *		COSB, NA, coef, cosA, sinA, sinD, cosD

C	BETAPRIME, ALPHA - SUN DIRECTION IN DEGREES
C	MASS - IN KG
C	ORB - IN KM, DEG, SEC
C	DT - IN SEC
C	AREA - IN M**2
C	CR - DIMENSIONLESS
C	C - SPEED OF LIGHT IN M/SEC

	CHARACTER*24 ATIME, AJUL4
	LOGICAL ERROR

	DATA AJUL4/'04-JUL-1995 00:00:00'/

	DOUBLE PRECISION PI
	PARAMETER (Pi = 3.14159265358979)

      double precision AE                       ! in kilomters
      double precision earth_freq               ! radians / second
      double precision earth_rate               ! meters / day
      double precision MU                       ! km**3/sec
      double precision GMMOON                   ! km**3/sec
      double precision GMSUN                    ! km**3/sec
      double precision sid_day                  ! seconds
      double precision deg_to_km                ! kilometers/deg
      double precision flat		
	double precision c

      common / physical_constants /
     &    AE,         earth_freq,  earth_rate,    MU,
     &    GMMOON,          GMSUN,     sid_day,    deg_to_km, flat, 
     &	  c

C
C	COMPUTE PHASE ANGLE & SOLAR RADIATION FLUX AT EARTH
C

	PHASE = 2.0*PI*( CH2JD(ATIME, .TRUE., ERROR) - 
     &			CH2JD(AJUL4, .TRUE., ERROR))/365.25
C
C	SRP IN WATTS/METER**2

	S = 1358/(1.0004+0.0334*DCOS(PHASE))
C
C	NET TOTAL FORCE IN NEWTONS

	F = S * AREA * CR / C

C
C	COMPUTE ANGULAR DURATION OF SHADOW, ASSUMING A CIRCULAR ORBIT
C	
	CMIN = DSQRT(1-(AE/ORB(1))**2)	

	COSB = DCOS(BETAPRIME*PI/180.0)
	cosA = dcos(alpha*pi/180.0)
        sinA = dsin(alpha*pi/180.0)

	IF (COSB.LT.CMIN) THEN
	   DELTA = PI
	ELSE
	   DELTA = PI - DACOS (CMIN/cosb)
	ENDIF

	sind = dsin(delta)
        cosd = dcos(delta)
C
C	compute velocity in meters/second
C
	na = SQRT(MU/ORB(1))*1000.0

	coef = F*cosb/(2*pi*mass*na)
C
C	compute solar radiation pressure rates on frozen orbit
C	
	dedt = coef*(-2.0*delta*sina+cosa*sina*sind*(4*cosa*cosd))
	dwdt = -(coef/orb(2))*cosa*(delta - 3.0*sind*cosd)

C	write(8,*) atime, ' PHASE = ', phase, ' S = ', S, 
C     *		' area = ', area, ' mass = ', mass, ' F = ', 
C     *		F, ' c = ', c, ' beta = ', betaprime
C     *		, ' alpha = ', alpha, ' delta = ', 
C     *		delta, ' dedt = ', dedt, ' dwdt = ', dwdt, 
C     *		' na = ', na, ' dt = ', dt
C
C	write(8,*) atime,' ',dedt*dt,' ', dwdt*dt*180.0/Pi
C	if (orb(2).ne.1) stop

C
C	update orbit
C
	orb(2) = orb(2) + dedt*dt
	orb(5) = orb(5) + dwdt*(180.0/Pi)*dt

	RETURN
	END