      SUBROUTINE FINDGB
C
C      SUBROUTINE FINDGB CALCULATES THE GRAVITY GRADIENT MATRIX,
C             SUNLINE VECTOR FOR INTERNAL ORBIT GENERATOR, AND
C             THE SOLAR ILLUMINATION FACTOR
C
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 J2,J3,J4,J22
C
C
      COMMON/CNBODY/ J2,J3,J4,J22,ZJ20,ZMU,WWO,FLAT,AEARTH
C
      COMMON/CONSTS/ PI,TWOPI,RADIAN
C
      COMMON/CSTVAL/ TSTART,ZL0(10),ZL1(10),ZL2(10),ZLA(10)
C
      COMMON/IMAIN1/ IDATE,LSAVE,INOPT,IPLOT,NUMEQS,IPLTPE,IORB,ITAPE
C
      COMMON/IPOOL1/ IGRAV,IDAMP,IK,K1,ITIM,IAB,IAPS,IBB,IBPS,NK(10),
     .               LK(10),LLK(10)
C
      COMMON/ISHADE/ IPLANS,ISATSH,IWRTTF
C
      COMMON/OUTTWO/ SOLILL,EPSERR
C
      COMMON/PRESUR/ DTOO,POO
C
      COMMON/RPOOL1/ RHOK(10),TIME,SA(3,3),FM1(3,3),ZLK(10),OMEG(3),
     .               ZLKP(10),ZLKDP(10),CMAT(3,3),GBAR(3,3),YBCM(3),
     .               ZBZK(3,10),FCM(3,3),DTO,PHID,PHI
C
      COMMON/RPOOL2/ PO,SD(3),DAN(3,10),DBN(3,10),CFMT(3,3),DIY1(3),
     .               SD1(3),DT1,P1,AERO,DTO1,YIZK(3),PO1
C
      COMMON /RSUNCL/ WE,TVER,ECLPTC
C
      COMMON/VECTRS/ XSAT(3),XSATDT(3),AD(3)
C
      COMMON/ORBNEW/ROD(3),VOD(3),TST
C
      DIMENSION XSAT1(3),XSAT2(3)
C
C
      DATA CONV/7.905355D0/
      IF(IORB.GT.0) GO TO 5
      IF(TIME-TSTART.GT.0) GO TO 4
C
      DO 3 I=1,3
      ROD(I)=XSAT(I)
    3 VOD(I)=XSATDT(I)
      TST=TSTART
    4 CALL XFIND(ROD,VOD,ZMU,TST   ,TIME,XSAT,XSATDT)
      GO TO 559
 5    CONTINUE
      TORB=TIME
      CALL TCNVRT(DDATE,MM,SEC,TLAST,TORB,2)
      CALL DTAPRE(ITAPE,IERR,DDATE,SEC,RX,RY,RZ,VX,VY,VZ)
      XSAT(1)=RX*AEARTH
      XSAT(2)=RY*AEARTH
      XSAT(3)=RZ*AEARTH
      XSATDT(1)=VX*CONV
      XSATDT(2)=VY*CONV
      XSATDT(3)=VZ*CONV
      IF(IERR.LT.3) GO TO 2000
      IF(IERR.LT.6) GO TO 3000
 559  CONTINUE
      RADSQR=0
      DO 10 I=1,3
   10 RADSQR=RADSQR + XSAT(I)*XSAT(I)
      RADIUS=DSQRT(RADSQR)
C
      DO 20 I=1,3
      XSAT1(I)=XSAT(I)/RADIUS
   20 XSAT2(I)=XSAT1(I)*XSAT1(I)
C
      DO 30 I=1,3
      DO 30 J=1,3
   30 GBAR(I,J)=0.D0
C
      IF(IGRAV.EQ.0) GO TO 60
C
      RAD3RD=RADIUS*RADSQR
      RAD4TH=RADIUS*RAD3RD
      RAD5TH=RADIUS*RAD4TH
C
      ZJ2OMU=-(ZJ20*ZMU*AEARTH*AEARTH)/(2.D0*RAD5TH)
      ZMUR3=-ZMU/RAD3RD
C
      DO 40 I=1,2
C
      GBAR(I,I)=  ZMUR3*(1.D0 - 3.D0*XSAT2(I))
     .          + ZJ2OMU*(3.D0 - 15.D0*(XSAT2(I) + XSAT2(3))
     .                  + 105.D0*XSAT2(I)*XSAT2(3))
C
   40 GBAR(3,I)=- 3.D0*ZMUR3*XSAT1(I)*XSAT1(3)
     .          + ZJ2OMU*(105.D0*XSAT1(I)*XSAT1(3)*XSAT2(3)
     .                  - 45.D0*XSAT1(I)*XSAT1(3))
C
C
C
      X21=XSAT1(1)*XSAT1(2)
      GBAR(2,1)=- 3.D0*ZMUR3*X21 + ZJ2OMU*(105.D0*X21*XSAT2(3)
     .          - 15.D0*X21)
C
      X34=XSAT2(3)*XSAT2(3)
C
      GBAR(3,3)=  ZMUR3*(1.D0 - 3.D0*XSAT2(3))
     .          + ZJ2OMU*(9.D0 - 90.D0*XSAT2(3) + 105.D0*X34)
C
C
      DO 50 I=1,3
      DO 50 J=I,3
   50 GBAR(I,J)=GBAR(J,I)
C
C                       SOLAR ILLUMINATION
C
   60 COSC=0.D0
C
      EPS=ECLPTC*RADIAN
      CEPS=DCOS(EPS)
      SEPS=DSIN(EPS)
      WS=WE*(TIME-TVER)*RADIAN
      SINWS=DSIN(WS)
      COSWS=DCOS(WS)
      SD(1)=-COSWS
      SD(2)=-SINWS*CEPS
      SD(3)=-SINWS*SEPS
C
   65 DO 70 I=1,3
   70 COSC=COSC + SD(I)*XSAT1(I)
C
      IF(COSC.LE.0.D0) GO TO 80
C
      SINC=DSQRT(1.D0 - COSC*COSC)
      RADCOS=RADIUS*COSC
      RADSIN=RADIUS*SINC - AEARTH
C
      H1=RADSIN + 0.00461D0*RADCOS
C
      IF(H1.LE.0.D0) GO TO 90
C
      H2=RADSIN - 0.004695D0*RADCOS
C
      IF(H2.GE.0.D0) GO TO 80
C
      H3=0.004653D0*RADCOS
      H4=1.D0 + H2/H3
C
      SOLILL=1.D0 - DARCOS(H4)/PI + H4/PI*DSQRT(1.D0 - H4*H4)
      GO TO 100
   80 SOLILL=1.D0
      GO TO 100
   90 SOLILL=0.D0
  100 DTO=SOLILL*DTOO
      PO=SOLILL*POO
      DT1=0.D0
C
      IF(IPLANS.EQ.0.AND.ISATSH.EQ.0) RETURN
      CALL SHADEP(SOLILL)
C
      RETURN
 2000 WRITE(6,2001)
 2001 FORMAT(' I/O ERROR ON EPHERMIS FILE')
 3000 WRITE(6,3001)
 3001 FORMAT(' INFORMATION NOT AVAILABLE FOR REQUESTED TIME')
      RETURN
      END
