subroutine plotmod C C Plots faked visibilities. Adapted from SCHEDULE for the binafit software C C C. A. Hummel 15 September 1991 C include 'constnts.inc' include 'data.inc' include 'model.inc' C integer*4 maxbase,niut parameter(maxbase=38) parameter(niut=48) real*8 modvis2(niut,mxfilt),modphas(niut,mxfilt),mhours(niut) real*8 lambda(niut,mxfilt),dlambda(niut,mxfilt) real*8 u(niut,mxfilt),v(niut,mxfilt),seprad,parad real*4 y(niut),t(niut) real*4 timmax,timmin integer*4 idate,mdate(niut),mbaseline(niut),mstars(niut) logical visible,doamp,dophase character header*80,device*4,plotfile*15 integer*4 bnum,spectrum,iut,firstiut,lastiut real *4 bx,by,bz,bcon,b(maxbase,4) REAL *4 zd,zlim,ha(niut),ra,dec real *4 delay,dlimit REAL *4 azi(niut),alt(niut),sinalt,sinazi,cosalt,cosazi real *4 st(niut),ut(niut),jddiff real *4 rasun,decsun,altsun,hasun,zsunlim REAL *8 tjd,gmst,jd,ram,decm REAL *4 pmra,pmdec,par,rv,mag,bv,angdia CHARACTER *8 name write(outc,'(a,$)')' Amplitude (a) or phase (p) ? ' read(inc,2)device doamp=.false. dophase=.false. if(device.eq.'a')then doamp=.true. else dophase=.true. endif write(outc,*)'Screen (s) or hardcopy (h,p,v) or exit (n) ?' read(inc,2)device 2 format(a4) if(device.eq.'h'.or.device.eq.'H')then plotfile='plotmod.lpm'//'/lpm' iy=4 csize=2.5 else if(device.eq.'p'.or.device.eq.'P')then plotfile='plotmod.ps'//'/ps' iy=4 csize=2.5 else if(device.eq.'v'.or.device.eq.'V')then plotfile='plotmod.ps'//'/vps' iy=4 csize=2.5 else if(device.eq.'s'.or.device.eq.'S')then plotfile='/xs' iy=3 csize=2.5 else return endif c Get date for which fake data should be produced write(outc,'(a,$)')' Enter date (yyyy mm dd): ' read(inc,*)iyear,imonth,iday idate=(iyear-1900)*10000+imonth*100+iday call julda(iyear,imonth,iday,tjd,gmst) gmst=gmst/3600.D0 c Read baseline file to get coordinates OPEN(UNIT=52,FILE='/home/chummel/mark3/binafit/baseline.dat', $ IOSTAT=IERR ) IF ( IERR .NE. 0 ) THEN WRITE(outc,*) '** Cannot open baseline.dat ' return ENDIF READ(52,'(A/A)') data,data DO I = 1, MAXBASE READ (52,'(I5,4F15.6)') BNUM, BX, BY, BZ, BCON if(bnum.ne.i)then write(outc,*)' Baseline file corrupt ! Return now !' return endif b(i,1)=bx b(i,2)=by b(i,3)=bz b(i,4)=bcon ENDDO close(52) c Find star information in catalog CALL CATALOG ( STARNO, NAME, RAM, DECM, PMRA, PMDEC, $ PAR, RV, MAG, BV, SPECTRUM, ANGDIA ) ra=real(ram) dec=real(decm) c This is the position of the sun jddiff=real(jd(iday,imonth,iyear)-jd(21,3,iyear)) rasun=jddiff*24./365. decsun=sin(jddiff*2.*pi/365.)*23.3 c Now find out at which times the star is observable C (1) points toward south point on horizon C (2) points toward east point on horizon C (3) points toward zenith C C Z(1)=0 Z(2)=0 Z(3)=1 Z = zenith C C Zenith angle limit ZLIM = 60. c Zenith angle limit for the sun ZSUNLIM=100. firstiut=1 lastiut=niut visible=.false. do iut=firstiut,lastiut utstep=24./niut ut(iut)=0.+(iut-1)*utstep st(iut)=real(gmst)+ut(iut)-8. ! 8 = time difference PST to GMST C HA(iut) = 15.*(ST(iut) - RA) IF ( HA(iut) .LT. -180. ) HA(iut) = HA(iut) + 360. IF ( HA(iut) .GT. 180. ) HA(iut) = HA(iut) - 360. SINALT = SIN(LAT*DEGRAD)*SIN(DEC*DEGRAD)+ $ COS(LAT*DEGRAD)*COS(DEC*DEGRAD)*COS(HA(iut)*DEGRAD) COSALT = COS( ASIN ( SINALT ) ) ALT(iut)= ASIN( SINALT )/DEGRAD SINAZI = COS(DEC*DEGRAD)*SIN(HA(iut)*DEGRAD) $ /COS(ALT(iut)*DEGRAD) COSAZI =(SIN(DEC*DEGRAD)-SINALT * SIN(LAT*DEGRAD))/ $ (COSALT * COS(LAT*DEGRAD)) AZI(iut)= ATAN2 ( SINAZI, COSAZI ) / DEGRAD IF ( AZI(iut) .LT. 0. ) AZI(iut) = 360. + AZI(iut) C zenith distance ZD = 90. - ALT(iut) c Now check, whether the sun is below the horizon HASUN = 15.*(ST(iut) - RASUN) IF ( HASUN .LT. -180. ) HASUN = HASUN + 360. IF ( HASUN .GT. 180. ) HASUN = HASUN - 360. SINALT = SIN(LAT*DEGRAD)*SIN(DECSUN*DEGRAD)+ $ COS(LAT*DEGRAD)*COS(DECSUN*DEGRAD)*COS(HASUN*DEGRAD) ALTSUN = ASIN( SINALT )/DEGRAD c zenith distance of sun ZDSUN = 90. - ALTSUN IF (ZD.LT.ZLIM.and.zdsun.gt.zsunlim) THEN visible=.true. if(firstiut.eq.1)firstiut=iut lastiut=iut endif enddo c Now loop over baselines and UT's DLIMIT = 9.5 do ig=1,nfilt ! initialize filter information wavehi=mfilter(ig)+25./2. wavelo=mfilter(ig)-25./2. wavest=(wavehi-wavelo)/(nusteps-1) do l=1,nusteps ftcurve(ig,l,1)=dble(l-1)*dble(wavest)+dble(wavelo) ftcurve(ig,l,2)=100.d0 ! 100% Transmission enddo enddo call pgbegin(0,plotfile,3,iy) do 100 ibase=1,maxbase if(b(ibase,1).eq.0)goto 100 k=1 timmin=ut(firstiut) timmax=ut(lastiut) do 110 iut=firstiut,lastiut C SH = SIN( HA(iut)*DEGRAD ) CH = COS( HA(iut)*DEGRAD ) SD = SIN( DEC*DEGRAD ) CD = COS( DEC*DEGRAD ) DELAY= B(ibase,4) + B(ibase,3)*SD $ - B(ibase,2)*CD*SH + B(ibase,1)*CD*CH IF ( ABS(DELAY) .LT. DLIMIT .and. visible)then mstars(k)=starno mdate(k)=idate mhours(k)=dble(ut(iut)) timmin=min(timmin,ut(iut)) timmax=max(timmax,ut(iut)) mbaseline(k)=ibase do j=1,nfilt lambda(k,j)=mfilter(j) dlambda(k,j)=j enddo k=k+1 ENDIF 110 continue ! Look for next time c Plot the model data, if you have more than 4 hours k=k-1 if(utstep*k.gt.3)then call calcuv(k,nfilt,niut,mxfilt,mstars,mdate,mbaseline, & mhours,lambda,u,v) mode=0 call binavis(k,nfilt,niut,mxfilt,mdate,mbaseline,mhours, & lambda,dlambda,modvis2,modphas,u,v, & seprad,parad,mode,ftcurve) pmax=0. do iid=1,k do iig=1,nfilt modphas(iid,iig)=modphas(iid,iig)/degrad pmax=max(pmax,abs(real(modphas(iid,iig)))) enddo enddo call pgsch(csize) do ig=1,nfilt v2max=0. do iid=1,k v2max=max(v2max,real(modvis2(iid,ig))) enddo v2max=v2max*1.2 call pgsci(5) if(doamp) call pgenv(timmin,timmax,0.,v2max,0,1) if(dophase)call pgenv(timmin,timmax,-pmax,pmax,0,0) write(header,3) & starno,idate,ibase,int(lambda(1,ig)) 3 format('FK',i4,1x,i6,1x,' B.l. ',i2,1x,i3,' nm') csize=csize+1 call pgsch(csize) call pglabel('','',header) csize=csize-1 call pgsch(csize) if(ig.eq.1)call pgsci(2) if(ig.eq.2)call pgsci(3) if(ig.eq.3)call pgsci(11) if(doamp)then do n=1,k y(1)=(real(modvis2(n,ig))*1.1)+0.002 y(2)=(real(modvis2(n,ig))*0.9)-0.002 t(1)=real(mhours(n)) t(2)=real(mhours(n)) call pgline(2,t,y) enddo else do n=1,k y(n)=real(modphas(n,ig)) t(n)=real(mhours(n)) enddo call pgline(k,t,y) endif enddo ! End of ig-loop endif 100 continue ! Next baseline call pgend return end c*********************************************************************** SUBROUTINE CATALOG( ISTAR, NAME, RAM, DECM, PMRA, PMDEC, PAR, $ RV, MAG, BV, SPECTRUM, ANGDIA ) C C CATALOG finds a star in the FK5 catalog. Returns star name ( NAME) C position ( RAM, DECM ) in fractional hours and degrees C proper motion ( PMRA, PMDEC ) parallax (PAR), C radial velocity ( RV), V magnitude (MAG) C C written and tested, 7 Dec, 1986 David Mozurkewich C FK5 catalog added, 6 Mar, 1987 C IMPLICIT NONE SAVE INTEGER*4 ISTAR, LCAT, I, IERR, INAME(2), STARNO, SPECTRUM REAL*8 RAM, DECM REAL*4 PAR, RV, MAG, PMRA, PMDEC, RAD, ANGDIA, BV CHARACTER *8 NAME LOGICAL FIRST DATA FIRST / .TRUE. / DATA LCAT, RAD / 53, 57.29577951D0 / C----------------------------------------------------------------------- IF ( FIRST ) THEN OPEN ( UNIT = LCAT, FILE='/home/chummel/mark3/binafit/fk5.bin', $ RECL=60, $ FORM='UNFORMATTED', ACCESS='DIRECT', STATUS='OLD', $ IOSTAT=IERR ) IF ( IERR .NE. 0 ) THEN WRITE(0,*) ' ERROR', IERR, ' OPENING ', $ ' X:\\FK5.BIN ' STOP END IF FIRST = .FALSE. END IF C IF ( (ISTAR .LE. 0) .OR. (ISTAR.GT.2000) ) THEN WRITE(0,*) ' Bad FK5 number ', ISTAR RAM = 0.0 DECM = 0.0 ELSE READ( LCAT,REC=ISTAR) STARNO, INAME(1), INAME(2), $ RAM, DECM, PMRA, PMDEC, PAR, RV, MAG, BV, SPECTRUM, ANGDIA IF ( STARNO .NE. ISTAR ) THEN WRITE(0,*) ' CATALOG is corrupt! ' WRITE(0,*) ' Record ', ISTAR, ' contains star ', STARNO WRITE(0,*) ' CATALOG bombs ' STOP END IF WRITE(NAME,'(2A4)') INAME(1), INAME(2) END IF RETURN 200 CONTINUE WRITE(0,1021) I, IERR STOP 1021 FORMAT ( ' fatal error while reading FK5.bin catalog '/ $ ' error number:', I5, ' record number:', I5/ ' too bad !') END