subroutine plotorb(chisq) c c The orbit plotting routine for the binafit software c c C. A. Hummel 15 September 1991 include 'constnts.inc' include 'data.inc' include 'model.inc' real*8 chisq,jd,cjul,djul,tjdstep,eposc,equinox integer*4 maxd,maxw,maxie,mxpos parameter (maxd=360) ! number of model orbit points parameter (maxw=12) ! number of windows parameter (maxie=400) ! number of steps in ellipse line integral parameter (mxpos=200) ! number of binary positions in plot character prognm*8,versn*3,hist*80,starname*30 character text*7,starn*8,orbsense*10 character fknam1*1,fknam2*2,fknam3*3,fknam4*4 character plotfile*15,plotname*8,device*4,command*1 character*1 digits(13) character*3 months(12),congens(88),greek(24),greeks(24) character*20 congen(88) character*80 title(19),data real*4 ra4(maxd),dec4(maxd),hra4(4),hdec4(4) real*4 minra4,mindec4,maxra4,maxdec4,rarnge,decrnge,range real*4 minra4o,mindec4o,maxra4o,maxdec4o real*4 wminra4(maxw),wmaxra4(maxw),wmindec4(maxw),wmaxdec4(maxw) integer*4 orbsymb logical plotcurr,query,zoom,ploterr,multiw,nodata,ask logical hardcopy,repeat,foundnum,dotitle,dowindow,doparbox c These are for plotting error ellipses real*4 dx(mxepl),dy(mxepl) real*4 ep,er,emaj,emin,epa,esqr,phi,radius,phistep,lell,dlell c These are for command 'i' real*4 raorb(mxpos),decorb(mxpos),chiorb(mxpos) integer*4 dateorb(mxpos) character*6 orbdate,orbchi c This is for the parameter box real*4 hparm(12),hparmerr(12) c These are for writing bintable.tex character*3 months2(12) real*4 blength(38) c These are for subroutine true2app real*8 a(7),dyda(7,2),ra,dec,seprad,parad,nparad integer*4 pr pr=6 a(1)=a_true*mas2rad a(2)=e_true a(3)=period a(4)=epoch a(5)=inclin*degrad a(6)=asc_node*degrad a(7)=arg_peri*degrad data blength/12.0,12.0, 3.0, 4.3, 5.3, 6.6, 6.9, 7.5, 8.2, 9.8, & 10.6,11.4,10.0, 0.0, 0.0, 0.0, 0.0,15.2, 0.0,15.4, & 0.0, 0.0, 0.0, 0.0,19.1,19.4, 0.0,19.7, 0.0, 0.0, & 23.1,23.3,23.6,25.2,27.0,27.6,29.1,31.5/ data months2/'Jan','Feb','Mar','Apr','May','Jun', & 'Jul','Aug','Sep','Oct','Nov','Dec'/ data months/'JAN','FEB','MAR','APR','MAY','JUN', & 'JUL','AUG','SEP','OCT','NOV','DEC'/ data congen/'Andromedae','Antliae','Apodis','Aquarii', &'Aquilae','Arae','Arietis','Aurigae','Bootis','Caeli', &'Camelopardi','Cancri','Canum Venaticorum','Canis Majoris', &'Canis Minoris','Capricorni','Carinae','Cassiopeiae','Centauri', &'Cephei','Ceti','Chamaeleontis','Circini','Columbae', &'Comae Berenices','Coronae Austrinae','Coronae Borealis', &'Corvi','Crateris','Crucis','Cygni','Delphini','Doradus', &'Draconis','Equulei','Eridani','Fornacis','Geminorum','Gruis', &'Herculis','Horologii','Hydrae','Hydri','Indi','Lacertae', &'Leonis','Leonis Minoris','Leporis','Librae','Lupi','Lyncis', &'Lyrae','Mensae','Microscopii','Monocerotis','Muscae','Normae', &'Octantis','Ophiuchi','Orionis','Pavonis','Pegasi','Persei', &'Phoenicis','Pictoris','Piscium','Piscis Austrini','Puppis', &'Pyxidis','Reticuli','Sagittae','Sagittarii','Scorpii', &'Sculptoris','Scuti','Serpentis','Sextantis','Tauri', &'Telescopii','Trianguli','Trianguli Australis','Tucanae', &'Ursae Majoris','Ursae Minoris','Velorum','Virginis', &'Volantis','Vulpeculae'/ c data congens/'AND','ANT','APS','AQR','AQL','ARA','ARI','AUR', 1'BOO','CAE','CAM','CNC','CVN','CMA','CMI','CAP','CAR','CAS', 2'CEN','CEP','CET','CHA','CIR','COL','COM','CRA','CRB','CRV', 3'CRT','CRU','CYG','DEL','DOR','DRA','EQU','ERI','FOR','GEM', 4'GRU','HER','HOR','HYA','HYI','IND','LAC','LEO','LMI','LEP', 5'LIB','LUP','LYN','LYR','MEN','MIC','MON','MUS','NOR','OCT', 6'OPH','ORI','PAV','PEG','PER','PHE','PIC','PSC','PSA','PUP', 7'PYX','RET','SGE','SGR','SCO','SCL','SCT','SER','SEX','TAU', 8'TEL','TRI','TRA','TUC','UMA','UMI','VEL','VIR','VOL','VUL'/ data greek/'\ga','\gb','\gg','\gd','\ge','\gz','\gy','\gh', 1'\gi','\gk','\gl','\gm','\gn','\gc','\go','\gp','\gr','\gs', 2'\gt','\gu','\gf','\gx','\gq','\gw'/ data greeks/'ALP','BET','GAM','DEL','EPS','ZET','ETA','THE', 1'IOT','KAP','LAM','MU ','NU ','XI ','OMI','PI ','RHO','SIG', 2'TAU','UPS','PHI','CHI','PSI','OME'/ data digits/'0','1','2','3','4','5','6','7','8','9','-','+','.'/ c c Initialize some variables c nw=0 fds=1.0 zero=0.1 orbsymb=1 hardcopy=.true. repeat=.false. dotitle=.true. doparbox=.false. ploterr=.true. plotcurr=.false. zoom=.false. ask=.false. dowindow=.false. if(ndata.eq.0)then nodata=.true. write(star,'(i5)')starno else nodata=.false. endif c write(outc,*)' Hit RETURN to continue, -n- to exit !' c read(inc,2)device device='s' 2 format(a4) 5 if(device.eq.'H')then write(outc,'(a,$)')' Name of plot file: ' read(inc,'(a)') plotname plotfile=plotname//'/lpm' else if(device.eq.'c')then plotfile='plotorb.cps'//'/cps' else if(device.eq.'P')then write(outc,'(a,$)')' Name of plot file: ' read(inc,'(a)') plotname plotfile=plotname//'/ps' else if(device.eq.'p')then plotfile='plotorb.ps'//'/ps' else if(device.eq.'V')then write(outc,'(a,$)')' Name of plot file: ' read(inc,'(a)') plotname plotfile=plotname//'/vps' else if(device.eq.'v')then plotfile='plotorb.ps'//'/vps' else if(device.eq.'n'.or.device.eq.'N')then return else plotfile='/xs' hardcopy=.false. endif c c The multi-windows option c iw=1 if(hardcopy.or.repeat)goto 6 if(query('Multi-windows ? (y/n) '))then multiw=.true. zoom=.true. write(outc,*)'Enter number of windows! (0=read from windows.dat)' read(inc,*)nw if(nw.eq.0)then open(9,file='windows.dat',status='old',err=12) ierr=0 do while(nw.lt.maxw.and.ierr.eq.0) nw=nw+1 read(9,*,iostat=ierr) & wmaxra4(nw),wminra4(nw),wmaxdec4(nw),wmindec4(nw) enddo if(ierr.ne.0)nw=nw-1 write(outc,*)'Read ',nw,' windows.' close(9) goto 13 endif 12 if(nw.eq.0)then write(outc,*)'*** No windows read. Return now!' return endif do i=1,nw write(outc,'(a,i1,a,$)')' Enter window ',i,' (+l -r +t -b):' read(inc,*)wmaxra4(i),wminra4(i),wmaxdec4(i),wmindec4(i) enddo 13 write(outc,*)'Enter number of windows in x and y dir.!' read(inc,*)nwx,nwy else multiw=.false. endif 6 continue c c Write nice greek star name c starname(1:5)=star(1:5) do i=1,24 j=index(star(1:3),greeks(i)) if(j.ne.0)then starname(1:3)=greek(i) goto 10 endif enddo 10 do i=1,88 j=index(star(6:8),congens(i)) if(j.ne.0)then starname(6:30)=congen(i) goto 20 endif enddo starname=star 20 continue c c Evaluate current date c if(.not.plotcurr)then call mkhist(prognm,versn,hist) do i=1,12 j=index(hist,months2(i)) if(j.ne.0)then read(hist(j-3:j-2),*)icday read(hist(j+4:j+7),*)icyear read(hist(j+9:j+10),*)ichour c add the difference of EST to UT ichour=ichour+5 if(ichour.ge.24)then ichour=ichour-24 icday=icday+1 endif icmonth=i goto 31 endif enddo endif 31 continue c cjul=(jd(icday,icmonth,icyear)-2440000.D0)+dble(ichour)/24.D0 text=' ' if(plotcurr)then if(iyear.lt.2000)then write(text(1:6),'(i2.2,i2.2,i2.2)')icyear-1900,icmonth,icday else write(text(1:6),'(i2.2,i2.2,i2.2)')icyear-2000,icmonth,icday endif endif c c Now calculate ellipse c tjdstep=period/dble(maxd) maxra4=0 minra4=0 maxdec4=0 mindec4=0 do i=1,maxd djul=dble(i-1)*tjdstep+cjul call true2app(djul,a,seprad,parad,ra,dec,dyda,pr) ra4(i)=real(ra/mas2rad) dec4(i)=real(dec/mas2rad) maxra4=max(maxra4,ra4(i)) minra4=min(minra4,ra4(i)) maxdec4=max(maxdec4,dec4(i)) mindec4=min(mindec4,dec4(i)) enddo c determine sense of orbital motion djul=djul+tjdstep call true2app(djul,a,seprad,nparad,ra,dec,dyda,pr) if(nparad.gt.parad)then sense=-1 else sense=1 endif c c Write model.dat which contains model r,theta-values c c Open MODEL.DAT file for output of r-theta values open(8,file='model.dat') if(inclin.gt.90.D0)then orbsense='RETROGRADE' else orbsense='PROGRADE' endif starn=star do i=1,8 if(starn(i:i).eq.' ')starn(i:i)='_' enddo write(8,3)starn,starno,period,orbsense 3 format(a8,1x,i4,2x,f8.2,2x,a10) c Now write MODEL.DAT file in ELLIPSE-format idate=0 eposc=0.D0 do id=1,ndata if(date(id).ne.date(id+1).or.baseline(id+1).ne.baseline(id) & .or.id+1.gt.ndata)then IDAY = MOD(DATE(ID), 100) IMONTH= MOD(DATE(ID)/100, 100) IYEAR = MOD(DATE(ID)/10000, 100) + 1900 ihour=8 ! Midnight on Mt Wilson djul=(jd(iday,imonth,iyear)-2440000.D0)+dble(ihour)/24.D0 call true2app(djul,a,seprad,parad,ra,dec,dyda,pr) write(8,4)seprad/mas2rad,parad/degrad, & zero,zero,zero,date(id),ihour*3600. 4 format(f6.2,1x,f7.2,2x,f6.3,1x,f6.3,1x,f6.1,2x,i6,1x,f6.0) idate=idate+1 eposc=eposc+djul endif enddo close(8) if(.not.nodata)then eposc=eposc/dble(idate)+2440000.D0 equinox=1950.D0+(eposc-2433282.5D0)/365.25D0 endif c c Now plot the box... c 11 continue !---------- continue plotting multi-windows here range=maxra4-minra4 maxra4=maxra4+range/10. minra4=minra4-range/10. range=maxdec4-mindec4 maxdec4=maxdec4+range/10. mindec4=mindec4-range/10. if(multiw)then maxra4=wmaxra4(iw) minra4=wminra4(iw) maxdec4=wmaxdec4(iw) mindec4=wmindec4(iw) elseif(dowindow)then maxra4=maxra4o minra4=minra4o maxdec4=maxdec4o mindec4=mindec4o endif maxra4o=maxra4 minra4o=minra4 maxdec4o=maxdec4 mindec4o=mindec4 rarnge=maxra4-minra4 ! the following is for constant decrnge=maxdec4-mindec4 ! shift calculation in any scale range=max(decrnge,rarnge) if(decrnge.eq.range)dsf=12. if(rarnge .eq.range)dsf=19. if(rarnge/decrnge.lt.1.5.and.dsf.eq.19.)then dsf=12. range=decrnge endif ddshift=fds*0.2*range/dsf if(multiw.and.iw.eq.1)then call pgbegin(0,plotfile,nwx,nwy) elseif(.not.multiw.and..not.repeat)then call pgbegin(0,plotfile,1,1) call pgask(ask) endif CALL PGSCI(5) ! 5 = CYAN call pgslw(1) if(multiw)then call pgsch(1.0) call pgadvance call pgvstand call pgwnad(maxra4,minra4,mindec4,maxdec4) call pgsch(2.5) call pgbox('BCNST',0.0,0,'BCNST',0.0,0) else call pgsch(1.0) CALL PGENV(maxra4,minra4,mindec4,maxdec4,1,1) endif C CALL EGA_RESTORE_DEFAULT CALL PGSCI(1) ! 1 = WHITE C C UPDATE AND WRITE TITLES C c ....Clear titles do i=1,19 title(i)=' '// $ ' ' enddo c ....Label axes and print header write(title(1),1000)starname 1000 format(' Apparent orbit of ',a30) if(.not.multiw)then !--------------- no titles with multi-windows call pgsci(7) ! 7 = yellow if(dotitle)then call pgsch(1.5) CALL PGMTEXT('T',1.6,0.5,0.5,TITLE(1)) endif call pgsci(5) call pgsch(1.0) title(2)='R.A. SEPARATION (mas)' title(3)='DEC. SEPARATION (mas)' CALL PGMTEXT('B',2.3,0.5,0.5,TITLE(2)) CALL PGMTEXT('L',2.3,0.5,0.5,TITLE(3)) call pgsci(1) c ....Write parameter box frame write(title(4),*) ' a \(0726) +/- mas' write(title(5),*) ' e \(0726) +/- ' write(title(6),*) ' i \(0726) +/- deg' write(title(7),*) ' \gW \(0726) +/- deg' write(title(8),*) ' P \(0726) +/- d ' write(title(9),*) ' \gw \(0726) +/- deg' write(title(10),*)' T \(0726) +/- JD ' write(title(11),*)' D1 \(0726) +/- mas' write(title(12),*)' D2 \(0726) +/- mas' write(title(13),1001)int(mfilter(1)) write(title(14),1001)int(mfilter(2)) write(title(15),1001)int(mfilter(3)) 1001 format('\gDm\d',i3,'nm\u\(0726) +/- mag') c c Now get parameter uncertainties, since we want to include c them in the parameter list and to use them for the number of c significant digits. c hparm(1)=real(a_true) hparmerr(1)=real(parmerr(10)) hparm(2)=real(e_true) hparmerr(2)=real(parmerr(11)) hparm(3)=real(inclin) hparmerr(3)=real(parmerr(14)) hparm(4)=real(asc_node) hparmerr(4)=real(parmerr(15)) hparm(5)=real(period) hparmerr(5)=real(parmerr(12)) hparm(6)=real(arg_peri) hparmerr(6)=real(parmerr(16)) hparm(7)=real(epoch) hparmerr(7)=real(parmerr(13)) hparm(8)=real(diam(1,1)) hparmerr(8)=real(parmerr(1)) hparm(9)=real(diam(2,1)) hparmerr(9)=real(parmerr(2)) hparm(10)=real(delmag(2,1)) hparmerr(10)=real(parmerr(7)) hparm(11)=real(delmag(2,2)) hparmerr(11)=real(parmerr(8)) hparm(12)=real(delmag(2,3)) hparmerr(12)=real(parmerr(9)) c ....Determine number of significant digits (nsdig) and fill in c ....the rest of the parameter box do 50 iparm=1,12 if(iparm.ge.1.and.iparm.le.9)then ipos1=16 ipos2=23 ipos3=28 ipos4=35 endif if(iparm.gt.9)then ipos1=22 ipos2=29 ipos3=34 ipos4=41 endif if(iparm.eq.4.or.iparm.eq.6)then ipos1=17 ipos2=24 ipos3=29 ipos4=36 endif if(hparmerr(iparm).eq.0)goto 1113 c nsdig=nint(log10(abs(hparm(iparm))/hparmerr(iparm)) c & -log10(abs(hparm(iparm)))+0.5) nsdig=nint(log10(1./hparmerr(iparm))+0.5) if(nsdig.lt.0)nsdig=0 if(nsdig.gt.7)nsdig=7 goto(1110,1111,1112,1113,1114,1115,1116,1117)nsdig+1 1110 write(title(iparm+3)(ipos1:ipos2),1100)hparm(iparm) write(title(iparm+3)(ipos3:ipos4),1120)hparmerr(iparm) goto 50 1111 write(title(iparm+3)(ipos1:ipos2),1101)hparm(iparm) write(title(iparm+3)(ipos3:ipos4),1121)hparmerr(iparm) goto 50 1112 write(title(iparm+3)(ipos1:ipos2),1102)hparm(iparm) write(title(iparm+3)(ipos3:ipos4),1122)hparmerr(iparm) goto 50 1113 write(title(iparm+3)(ipos1:ipos2),1103)hparm(iparm) write(title(iparm+3)(ipos3:ipos4),1123)hparmerr(iparm) goto 50 1114 write(title(iparm+3)(ipos1:ipos2),1104)hparm(iparm) write(title(iparm+3)(ipos3:ipos4),1124)hparmerr(iparm) goto 50 1115 write(title(iparm+3)(ipos1:ipos2),1105)hparm(iparm) write(title(iparm+3)(ipos3:ipos4),1125)hparmerr(iparm) goto 50 1116 write(title(iparm+3)(ipos1:ipos2),1106)hparm(iparm) write(title(iparm+3)(ipos3:ipos4),1126)hparmerr(iparm) goto 50 1117 write(title(iparm+3)(ipos1:ipos2),1107)hparm(iparm) write(title(iparm+3)(ipos3:ipos4),1127)hparmerr(iparm) goto 50 1100 format(f8.0) 1101 format(f8.1) 1102 format(f8.2) 1103 format(f8.3) 1104 format(f8.4) 1105 format(f8.5) 1106 format(f8.6) 1107 format(f8.7) 1120 format(f6.0) 1121 format(f6.1) 1122 format(f6.2) 1123 format(f6.3) 1124 format(f6.4) 1125 format(f6.5) 1126 format(f6.5) 1127 format(f6.5) 50 continue c ....Write some additional information into the parameter box if(.not.nodata)then write(title(16),1002)real(chisq/ndata/nfilt) 1002 FORMAT(' \gx\u2\d/\gn = ',F5.2) write(title(17),1003)ndata*nfilt 1003 format(' # vis. = ',i4) write(title(18),1004)int(eposc) 1004 format(' Epoch of osculation: ',i7) write(title(19),1005)equinox 1005 format(' Mean Equinox: ',f6.1) else write(title(16),1002)real(chioc) endif endif !--------------- no titles with multi-windows -----------! call pgsch(1.5) C PLOT MODEL ORBIT POINTS ON THE SCREEN if(.not.zoom)then if(ploterr)then call pgsci(6) ! Magenta (red + blue) call pgpoint(maxd,ra4,dec4,orbsymb) else call pgsci(7) ! Yellow call pgpoint(maxd,ra4,dec4,orbsymb) endif else if(.not.multiw)then call pgpoint(maxd,ra4,dec4,orbsymb) endif call pgsci(6) call pgslw(1) call pgline(maxd,ra4,dec4) if(.not.multiw)then hra4(2)=ra4(1) hdec4(2)=dec4(1) hra4(1)=ra4(maxd) hdec4(1)=dec4(maxd) call pgline(2,hra4,hdec4) endif endif c c Plot r-theta positions of single nights, read from `FK5#'.dat c do j=1,mxpos ! initialize arrays for use with 'i' command raorb(j)=0. decorb(j)=0. dateorb(j)=0 enddo if(ploterr.or.nodata)then write(fknam4,'(i4)')starno if(log10(real(starno)).lt.1)then write(fknam1,'(i1)')starno open(4,file=fknam1//'.dat',status='old',err=42) else if(log10(real(starno)).lt.2)then write(fknam2,'(i2)')starno open(4,file=fknam2//'.dat',status='old',err=42) else if(log10(real(starno)).lt.3)then write(fknam3,'(i3)')starno open(4,file=fknam3//'.dat',status='old',err=42) else open(4,file=fknam4//'.dat',status='old',err=42) endif read(4,*,end=42) ! skip header isymb=1 c .....Determine, how many numbers are on a line in .DAT-file read(4,'(a80)',end=42)data foundnum=.false. innums=0 do 70 i=1,80 do j=1,13 if(data(i:i).eq.digits(j))then foundnum=.true. goto 70 endif enddo if(foundnum)then innums=innums+1 foundnum=.false. endif 70 continue rewind 4 read(4,*) ! skip header c......Prepare bintable output file open(9,file='bintable.tex') write(9,2001) 2001 format('\documentstyle[12pt,~cah/latex/sty/aasms]{article}') write(9,2002) 2002 format('\begin{document}') write(9,2003) 2003 format('\begin{table}') write(9,2004)starn 2004 format('\caption{Observation and result log for ',a30,'}') write(9,2005) 2005 format('\begin{tabular}{lcrrrrrrrr} \tableline \tableline') write(9,2006) 2006 format('&Bess. Yr.&&No.&$\rho$&$\theta$&$\sigma_{\rm maj}$&', & '$\sigma_{\rm min}$&$\varphi$&O-C\\ ') write(9,2007) 2007 format('UT Date&(-1900)&BL [m]&scans&[mas]&[deg]&[mas]&', & '[mas]&[deg]&[mas]\\ ') write(9,2008) 2008 format('(1)&(2)&(3)&(4)&(5)&(6)&(7)&(8)&(9)&(10)\\ \tableline') ityear=0 c chioc=0.0 c do j=1,mxpos if(innums.eq.9)then read(4,*,end=43)er,ep,emaj,emin,epa,idate,utsec,exc,eyc elseif(innums.eq.10)then read(4,*,end=43)er,ep,emaj,emin,epa,idate,utsec,exc,eyc,isymb endif if((j.eq.2.and.emin.gt.100000).or. $ (innums.ne.9.and.innums.ne.10))then write(outc,'(a,a,a)')' Sorry! Your format of ',fknam4//'.dat', $ ' is wrong!' write(outc,*)'BINAFIT expects: r, theta, e-maj, e-min, e-pa,', $ ' date, utsec, e-xcenter, e-ycenter, isymb (optional)' write(outc,*)'(e-prefix means: error ellipse)' close(4) goto 44 endif epa=epa+90. IDAY = MOD(IDATE, 100) IMONTH= MOD(IDATE/100, 100) IYEAR = MOD(IDATE/10000, 100) + 1900 ihour = 8 ! midnight on Mt Wilson = ref. hour for r-theta pos. ihour = utsec/3600 djul=(jd(iday,imonth,iyear)-2440000.D0)+dble(ihour)/24.D0 call true2app(djul,a,seprad,parad,ra,dec,dyda,pr) hra4(2)=real(seprad*sin(parad)/mas2rad) hdec4(2)=real(seprad*cos(parad)/mas2rad) epa=epa*real(degrad) hra4(1)=er*sin(ep*real(degrad)) hdec4(1)=er*cos(ep*real(degrad)) raorb(j)=hra4(1) ! use these later for the 'i' command decorb(j)=hdec4(1) dateorb(j)=idate esqr=1.0-emin*emin/(emaj*emaj) lambda=(emaj-emin)/(emaj+emin) lell=real(pi)*(emaj+emin)*(64.-3.*lambda**4)/(64.-16.*lambda**2) dlell=ddshift/fds/4. if(lell/dlell.lt.20)dlell=lell/20. if(lell/dlell.gt.200)dlell=lell/200. phi=0. phistep=360.0*real(degrad)/(maxie-1) lell=0. kk=0 radiuso=sqrt(emin*emin/(1.0-esqr)) sigoc=0.0 oc=sqrt((hra4(2)-hra4(1))**2+(hdec4(2)-hdec4(1))**2) do k=1,maxie phi=phi+phistep radius=sqrt(emin*emin/(1.0-esqr*cos(phi)**2)) lell=lell+sqrt(radius**2+radiuso**2 & -2.*radius*radiuso*cos(phistep)) radiuso=radius if(lell.gt.dlell)then if(kk.ge.mxepl)goto 80 kk=kk+1 lell=mod(lell,dlell) dx(kk)=(-(radius*cos(phi+epa)+exc)+hra4(1)) dy(kk)=(radius*sin(phi+epa)+eyc+hdec4(1)) sigoc=max(sigoc,((hra4(1)-hra4(2))*(dx(kk)-hra4(1)) & +(hdec4(1)-hdec4(2))*(dy(kk)-hdec4(1))) / oc) endif enddo 80 continue chioc=chioc+(oc/sigoc)**2 chiorb(j)=(oc/sigoc)**2 ! use this later for the 'i' command call pgsci(11) call pgslw(1) call pgsls(1) if(isymb.eq.27)then isymb=-1 call pgsls(4) call pgpoint(kk,dx,dy,isymb) else call pgline(kk,dx,dy) dx(2)=dx(1) ! now bridge gap dy(2)=dy(1) dx(1)=dx(kk) dy(1)=dy(kk) call pgline(2,dx,dy) endif call pgsci(7) call pgpoint(1,hra4,hdec4,isymb) call pgsci(6) call pgline(2,hra4,hdec4) call pgsls(1) c.......Now write to file bintable.tex if(iyear.ne.ityear)then write(9,2009)iyear 2009 format(i4,'&&&&&&&&&\\ ') ityear=iyear endif bessy=real((djul+24979.68648D0)/365.242198781D0) nscans=0 do id=1,ndata if(date(id).eq.idate)then ibase=baseline(id) nscans=nscans+1 endif enddo write(9,2010)months2(imonth),iday,bessy,blength(ibase),nscans, & er,ep,emaj,emin,epa/real(degrad)-90.,oc 2010 format(a3,1x,i2,'&',f7.4,'&',f4.1,'&',i3,'&',f6.2,'&',f6.2, & '&',f5.2,'&',f5.2,'&',f5.1,'&',f4.2,'\\ ') c enddo write(outc,*)' *** More dates in .DAT file!' goto 43 42 continue write(outc,'(a,a)')' *** File does not exist: ',fknam4//'.dat' write(outc,'(a)')' *** or is empty.' close(4) ploterr=.false. goto 44 43 close(4) write(9,2011) 2011 format('\tableline \tableline') write(9,2012) 2012 format('\end{tabular} \end{table} \end{document}') close(9) if(j-1-7.gt.0)then chioc=chioc/(j-1-7) else chioc=0 endif c open(7,file='plotorb.out') c do i=1,j-1 c write(7,*)chiorb(i) c enddo c close(7) else c c Alternatively, plot phases of the model corresponding to the visibilities c 44 call pgsci(7) do id=1,ndata IDAY = MOD(DATE(ID), 100) IMONTH= MOD(DATE(ID)/100, 100) IYEAR = MOD(DATE(ID)/10000, 100) + 1900 djul=(jd(iday,imonth,iyear)-2440000.D0)+hours(id)/24.D0 call true2app(djul,a,seprad,parad,ra,dec,dyda,pr) hra4(1)=real(ra/mas2rad) hdec4(1)=real(dec/mas2rad) call pgpoint(1,hra4,hdec4,17) enddo endif c c Now plot the current position of the binary call pgsch(1.5) call pgsci(6) ! Magenta hra4(1)=0. hdec4(1)=0. call true2app(cjul,a,seprad,parad,ra,dec,dyda,pr) hra4(2)=real(ra/mas2rad) hdec4(2)=real(dec/mas2rad) if(plotcurr)then call pgpoint(2,hra4,hdec4,18) call pgsls(2) call pgline(2,hra4,hdec4) call pgsls(1) endif c ...and a little arrow indicating the sense of rotation djul=cjul+tjdstep call true2app(djul,a,seprad,parad,ra,dec,dyda,pr) hra4(3)=real(ra/mas2rad) hdec4(3)=real(dec/mas2rad) angle=atan2(hra4(3)-hra4(2),hdec4(3)-hdec4(2))/real(degrad) call pgsch(1.0) if(.not.multiw) $ call pgptext(hra4(2),hdec4(2),angle,1.05,text//'\(2262)') c Now plot the semimajor axis connecting the periastron point c and the centre (the focus [i.e. the point (0,0)] should lie also c on that line). If the eccentricity is zero, plot the line of nodes. call pgsch(1.0) call pgsci(8) djul=epoch if(e_true.lt.0.001d0)djul=djul-arg_peri/360.d0*period call true2app(djul,a,seprad,parad,ra,dec,dyda,pr) hra4(1)=real(ra/mas2rad) hdec4(1)=real(dec/mas2rad) if(multiw)then call pgsch(2.5) if (hra4(1).gt.minra4 .and.hra4(1) .lt.maxra4 $ .and.hdec4(1).lt.maxdec4.and.hdec4(1).gt.mindec4) $ call pgptext(ra4(1),dec4(1),0.0,0.0,'T') else if(e_true.lt.0.001d0)then call pgptext(hra4(1),hdec4(1),0.0,0.0,'T\d0') else call pgptext(hra4(1),hdec4(1),0.0,0.0,'T') endif endif djul=djul+period/2.D0 call true2app(djul,a,seprad,parad,ra,dec,dyda,pr) hra4(2)=real(ra/mas2rad) hdec4(2)=real(dec/mas2rad) call pgline(2,hra4,hdec4) c ... Plot circles indicating the size of the stars c hra4(1)=0. c hdec4(1)=0. c phistep=360.0*real(degrad)/(mxepl-1) c do k=1,ncomp c phi=0. c radius=real(diam(k,1))/2. c do kk=1,mxepl c phi=phi+phistep c dx(kk)=radius*cos(phi)+hra4(k) c dy(kk)=radius*sin(phi)+hdec4(k) c enddo c call pgsls(1) c call pgslw(1) c call pgline(mxepl,dx,dy) c enddo iw=iw+1 ! multi-windows option: increase window counter if(hardcopy.and.iw.gt.nw)then if(.not.doparbox)then call pgend return endif dshift=0. if(device.eq.'v'.or.device.eq.'V')ddshift=ddshift/1.5 call pgsch(fds*0.5) do i=4,19 call pgptext(cx,cy-dshift,0.,0.,title(i)) dshift=dshift+ddshift enddo call pgend return endif c c Now plot cursor and read command c repeat=.false. if(.not.multiw.or.iw.gt.nw)then 60 continue write(outc,*)'Click plot window and enter command, ? for help:' command=' ' call pgcurse(cx,cy,command) c rewind inc if(command.eq.'?')then write(outc,*)'w Write parameter box' write(outc,*)'b Remove parameter box' write(outc,*)'p/P PostScript plot' write(outc,*)'v/V Vert. PS plot' write(outc,*)'c/C Color PS plot' write(outc,*)'o Enter new model orbit symbol' write(outc,*)'t Title toggle switch' write(outc,*)'d Position at spec. date' write(outc,*)'i Identify error ellipse' write(outc,*)'e Error ellipses (toggle)' write(outc,*)'l Plot orbit line (toggle)' write(outc,*)'z Enter new window' goto 60 elseif(command.eq.'q')then call pgend return elseif(command.eq.'c'.or.command.eq.'p'.or.command.eq.'v')then hardcopy=.true. device=command goto 5 elseif(command.eq.'C'.or.command.eq.'P'.or.command.eq.'V')then hardcopy=.true. device=command goto 5 elseif(command.eq.'o')then ! Get symbol for plotting model orbit points write(outc,'(a,$)') & ' Symbol for orbit points (T=tiny, D=dot, C=circle, O=other): ' read(inc,'(a)') command if(command.eq.'t' .or. command.eq.'T') then orbsymb=-1 else if(command.eq.'d' .or. command.eq.'D') then orbsymb=1 else if(command.eq.'c' .or. command.eq.'C') then orbsymb=22 else write(outc,'(a,$)') ' Type PGPLOT symbol number: ' read(inc,*) orbsymb end if elseif(command.eq.'t')then dotitle=.not.dotitle elseif(command.eq.'d')then write(outc,'(a,$)') & ' Enter the current date (yyyy mm dd [8 UT]): ' read(inc,*,end=30,err=30)icyear,icmonth,icday ichour=8 ! 8 UT = midnight at Mt Wilson if(icyear.eq.0.or.icmonth.eq.0.or.icday.eq.0)goto 30 plotcurr=.true. repeat=.true. goto 5 c 30 rewind inc 30 continue plotcurr=.false. elseif(command.eq.'e')then ploterr=.not.ploterr elseif(command.eq.'z')then write(outc,'(a,$)')' Enter window (+l -r +t -b): ' read(inc,*,end=40)maxra4o,minra4o,maxdec4o,mindec4o dowindow=.true. c 40 rewind inc 40 continue elseif(command.eq.'l')then zoom=.not.zoom elseif(command.eq.'s')then write(outc,'(a,$)')' Rescale par-box by: ' read(inc,*)factor fds=fds*factor ddshift=ddshift*factor goto 60 elseif(command.eq.'b'.or.command.eq.'w')then if(command.eq.'b')then call pgsci(0) doparbox=.false. else call pgsci(1) doparbox=.true. endif dshift=0. call pgsch(fds*0.5) do i=4,19 call pgptext(cx,cy-dshift,0.,0.,title(i)) dshift=dshift+ddshift enddo goto 60 elseif(command.eq.'i'.or.command.eq.'I')then j=1 if(dateorb(1).ne.0)then dist=(cx-raorb(j))**2+(cy-decorb(j))**2 idate=1 do while(dateorb(j).ne.0) if((cx-raorb(j))**2+(cy-decorb(j))**2.lt.dist)then dist=(cx-raorb(j))**2+(cy-decorb(j))**2 idate=j endif j=j+1 enddo cx=raorb(idate) cy=decorb(idate) write(orbdate,'(i6.6)')dateorb(idate) write(outc,*)'Date is : '//orbdate write(orbchi,'(f4.1)')chiorb(idate) c write(outc,*)'/ Chi^2 is: '//orbchi endif goto 60 endif ! end of command interpreter repeat=.true. goto 5 else ! continue with multi-window plotting goto 11 endif RETURN END