SUBROUTINE MODVIS
C Calculate the visibility of a model for plotting with BINPLOT;
C returns visibilities squared as MVIS2
C     J.T. Armstrong      12 March 1991

      INCLUDE 'BINFIT.INC'

      INTEGER*4 MXPOINT

      PARAMETER(MXPOINT=100)

      REAL*8 TMIN, TMAX, DELT, UT1, TDTOFF, DJUL, GAST, HA
      REAL*8 PMRA, PMDEC, PAR, RV, RAA, DECA, RA, DEC, RAM, DECM
      REAL*8 UMOD(MXPOINT,MXFILT), VMOD(MXPOINT,MXFILT)
      REAL*8 DMOD(MXPOINT,MXFILT)
      REAL*8 U(3),V(3),W(3),B(3)
      REAL*8 BX(MXBAS),BY(MXBAS),BZ(MXBAS)
      REAL*8 MVIS2(MXPOINT,MXFILT), T(MXPOINT)
      REAL*8 UVDIST, UVDIST1, UVANGLE, NU, CPHASE, BESJ0ARG
      REAL*8 BESJ1ARG, VISPHA, BESSJ1, DVISMOD
      REAL*8 JD, DOT
      INTEGER*4 INU, IDAY, IMON, IYEAR, NBODY, IERR, INDAT, INCAT
      COMPLEX*16 MVIS, MVISTOT(MXPOINT,MXFILT),DVIS

      COMMON /MOD_VIS/ MVIS2,T,TMIN,TMAX

      PARAMETER(INCAT=11)
      PARAMETER(INDAT=12)

      real*8 seprado,parado
      real*8 arg,alpha,beta,ul(mxmod,mxfilt)
      common/limb/ul

      EXTERNAL JD, DOT

C Convert angles (position angle, separation, diameter) to radians.

D     WRITE(OUTC,*) ' Convert angular quantities to radians'
      DO IC = 1, NCOMP
          DO IG = 1, NFILT
C              WRITE(OUTC,'(A,I1,A,I1,A,3(F10.5,2X))') 
C     1           ' MODVIS: Sep, PA, Diam (',IC,',',IG,'): ',
C     2           PARM(IC,IG,1),PARM(IC,IG,2),PARM(IC,IG,3)
              SEPRAD (IC,IG) = PARM(IC,IG,1)*MAS2RAD
              PARAD  (IC,IG) = PARM(IC,IG,2)*DEGRAD
              DIAMRAD(IC,IG) = PARM(IC,IG,3)*MAS2RAD
          END DO
      END DO
C Calculate the times for desired visibilities
      DELT = (TMAX-TMIN)/(MXPOINT-1)
      T(1) = TMIN
      DO I = 2,MXPOINT
          T(I) = T(I-1)+DELT
C          IF(I.GT.95) THEN
C              WRITE(OUTC,'(A,I3,A,F7.3)') ' T(',I,') = ',T(I)
C          END IF
      END DO
C Get baseline and star info
      OPEN(UNIT=INDAT,FILE='e:\\catalogs\\fk5\\BASELINE.DAT',
     $      STATUS='OLD',IOSTAT=IERR)
      READ(INDAT,'(1(/))') ! Skip header lines
      DO I=1,MXBAS
          READ(INDAT,'(I5,3(F15.6))',IOSTAT=IERR)
     1                  J, BX(I), BY(I), BZ(I)
      END DO
      CLOSE(INDAT)
      CALL READFK5( STARNO, STARNAME, RAM, DECM, PMRA, PMDEC, PAR, RV)
C Calculate u and v for those times
C
      DO J = 1, MXPOINT
          IF (J.EQ.1) THEN
              IDAY  = MOD(DATE,       100)
              IMON  = MOD(DATE/100,   100)
              IYEAR = MOD(DATE/10000, 100) + 1900
C             WRITE(OUTC,1300) J, DATE, IDAY, IMON, IYEAR
 1300         FORMAT( ' SCAN = ', I5, ' DATE = ', I10, 3I5,/ )
              DJUL  = JD ( IDAY, IMON, IYEAR )
C             WRITE(OUTC,'(A,F20.5)') ' DJUL = ',DJUL
              B(1)  = BX(BASELINE(J))     ! Baseline components in meters
              B(2)  = BY(BASELINE(J))
              B(3)  = BZ(BASELINE(J))
C             WRITE(OUTC,'(A,3(1X,G12.5))') ' B(1,2,3) (m): ',
C    1                                       B(1),B(2),B(3)
          END IF
          UT1  = T(J) / 24.D0
C         WRITE(OUTC,'(A,I4,A)') ' Calculating TJD(',J,')'
          TJD(J)  = DJUL + UT1 + TDTOFF
          NBODY = 3
          CALL APSTAR (TJD(J), NBODY, RAM, DECM, PMRA, PMDEC, PAR,
     $                         RV, RAA, DECA)
          CALL SIDTIM ( DJUL, UT1, 1, GAST )
          CALL DIURN ( TJD, GAST, RAA, DECA, RA, DEC )
          HA = 15.D0 *( GAST - RA ) - LONG
          IF ( HA .LT. -180 ) HA = HA + 360.D0
C         Components of u, v, w axes in local coordinate system
          U(1) = -SIN(HA*DEGRAD)
          U(2) = -COS(HA*DEGRAD)
          U(3) =  0.D0
          V(1) =  SIN(DEC*DEGRAD) * COS(HA*DEGRAD)
          V(2) = -SIN(DEC*DEGRAD) * SIN(HA*DEGRAD)
          V(3) = -COS(DEC*DEGRAD)
          W(1) = -COS(DEC*DEGRAD) * COS(HA*DEGRAD)
          W(2) =  COS(DEC*DEGRAD) * SIN(HA*DEGRAD)
          W(3) = -SIN(DEC*DEGRAD)
          DO IG = 1, NFILT
              UMOD(J,IG) = DOT( B, U ) / (LAMBDA0(IG)*1.0E-9)
              VMOD(J,IG) = DOT( B, V ) / (LAMBDA0(IG)*1.0E-9)
              DMOD(J,IG) = DOT( B, W ) / (LAMBDA0(IG)*1.0E-9)
          END DO
      END DO
D     WRITE(OUTC,*) ' U, V coordinates calculated'

C Now calculate the model visibilities:
      DO ID = 1, MXPOINT
D         WRITE(OUTC,'(A,I3,A)') ' Calculating vis''s at ',ID,'th time'
          DO IG = 1, NFILT
              MVISTOT(ID,IG) = CMPLX(0.0,0.0)
              DO IC = 1, NCOMP
c save seprad, parad, since their new values due to orbital motion
c are used
      seprado=seprad(ic,ig)
      parado=parad(ic,ig)
      if(ic.eq.2)then
      seprad(ic,ig)=seprad(ic,ig)+(tjd(id)-mepoch3)*dsep
      parad(ic,ig)=parad(ic,ig)+(tjd(id)-mepoch3)*dposang
      endif

      alpha=1.D0-ul(ic,ig)
      beta=ul(ic,ig)
          
                  UVDIST = DSQRT(UMOD(ID,IG)**2 + VMOD(ID,IG)**2)
                  UVDIST1 = UVDIST                    ! Diff. should be very small
                  UVANGLE = DATAN2(VMOD(ID,IG),UMOD(ID,IG))
C                 --Do the frequency integration:
                  MVIS = CMPLX(0.0D0, 0.0D0)
                  NU = NULO(IG) - 0.5 * DNU(IG)
                  DO INU = 1, NUSTEPS
                      NU       = NU + DNU(IG)
                      CPHASE   = TWOPI * SEPRAD(IC,IG) * (NU/NU0(IG))
     1                     * ( UMOD(ID,IG) * DSIN(PARAD(IC,IG))
     1                           + VMOD(ID,IG) * DCOS(PARAD(IC,IG)) )
                      BESJ1ARG = PI*DIAMRAD(IC,IG)*UVDIST1*(NU/NU0(IG))
c                     DVISMOD  = 2.0D0 * BESSJ1(BESJ1ARG)/BESJ1ARG
                      arg=besj1arg
             dvismod = (alpha*bessj1(arg)/arg+beta*sqrt(pi/2.D0)*
     &                  sqrt(2.D0/(pi*arg))*(sin(arg)/arg-cos(arg))/
     &                  sqrt(arg*arg*arg))/(alpha/2.D0+beta/3.D0)
                      DVIS     = CMPLX(DVISMOD*DCOS(CPHASE), 
     1                                           DVISMOD*DSIN(CPHASE))
                      MVIS     = MVIS + DVIS
                  END DO
C                 WRITE(OUTC,*) ' Finish frequency integration'
                  MVIS = MVIS / NUSTEPS
                  MVISTOT(ID,IG) = MVISTOT(ID,IG) + FLUX(IC,IG)*MVIS
      seprad(ic,ig)=seprado
      parad(ic,ig)=parado
              END DO
              MVIS2(ID,IG) = DBLE(MVISTOT(ID,IG))**2 
     1                       + DIMAG(MVISTOT(ID,IG))**2
          END DO
      END DO

      RETURN
      END