SUBROUTINE BINAVIS(ndata,nfilt,maxnd,maxnf,date,baseline,hours, & lambda,dlambda,modvis2,modphas,u,v, & mseprad,mparad,mode,ftcurve) C C Calculate the visibilities of a binary made up of two components, C with the phase center at the center of brightness, and separation C and position angle changing with time according to specified orbit. C Alternatively, take position angle and separation from sub. call. C save C include 'constnts.inc' include 'model.inc' integer*4 ndata,nfilt,maxnd,maxnf integer*4 date(maxnd),baseline(maxnd) real*8 hours(maxnd),lambda(maxnd,maxnf),dlambda(maxnd,maxnf) real*8 modvis2(maxnd,maxnf),modphas(maxnd,maxnf) real*8 u(maxnd,maxnf),v(maxnd,maxnf),ftcurve(mmxfilt,50,2) real*8 tjd,jd,djul,unm,vnm,nfactor real*8 seprad,parad,seprado,parado,mseprad,mparad real*8 bseprad(mxmod),bparad(mxmod) real*8 diamrad,color,flux,inlambda,flambda logical first,query,doimlc external jd C VISIBILITIES: COMPLEX*16 VIS, VISTOT, DVIS C Declarations for BWVIS part C Linear limb darkening variables real*8 arg,ul,alpha,beta C Planck curve real*8 cfactor,temp(mxmod),pcurve C U-V PLANE QUANTITIES: REAL*8 UVDIST C PHASES: REAL*8 CPHASE C INTEGRANDS AND VISIBILITIES: REAL*8 BESSJ1, DVISMOD c These are for subroutine true2app real*8 a(7),dyda(7,2),ra,dec integer*4 pr c Miscellaneous for the IMLC code integer*4 maxlc parameter (maxlc=1000) integer*4 clc(maxlc) real*8 xlc(maxlc),ylc(maxlc),flc(maxlc) common /imlcdata1/clc,ilc common /imlcdata2/xlc,ylc,flc real*8 flctot(mxmod),hxlc(maxlc),hylc(maxlc),rangle,ulam,vlam real*4 lcparms(54),lcoparms(54),spotparms(4,2,100) common /dspots/spotparms c data first/.true./ data doimlc/.false./ data lcoparms/54*0.0/ c c Initialize values for true2app 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 C if(ndata.lt.0.or.(first.and.starno.eq.193))then if(query('Use interacting binary code IMLC ? (y/n) '))then doimlc=.true. else doimlc=.false. endif endif if(doimlc)goto 10 C Save input values for seprad and parad seprado=mseprad parado=mparad if(starno.eq.193)then temp(1)=5700.d0 ! values for Capella Ab temp(2)=5000.d0 ! Aa elseif(starno.eq.227)then temp(1)=9000.d0 ! values for Bet Aur temp(2)=9000.d0 else temp(1)=6000.d0 temp(2)=temp(1) endif C----------------------------------------------------------------------- DO ID=1,NDATA C----------------------------------------------------------------------- C Set separation and position angle of binary C IF(DATE(ID).NE.DATE(ID-1).OR.ID.EQ.1)THEN IDAY = MOD(DATE(ID), 100) IMON = MOD(DATE(ID)/100, 100) IYEAR = MOD(DATE(ID)/10000, 100) + 1900 TJD = JD ( IDAY, IMON, IYEAR ) - 2440000.D0 ENDIF DJUL = TJD + HOURS(ID)/24.D0 call true2app(djul,a,seprad,parad,ra,dec,dyda,pr) if(mode.eq.1)then mseprad=mseprad+seprad mparad=mparad+parad ihour=8 ! Midnight on Mt Wilson djul=tjd+dble(ihour)/24.D0 call true2app(djul,a,seprad,parad,ra,dec,dyda,pr) mseprad=mseprad-seprad mparad =mparad-parad seprad =mseprad parad =mparad endif bseprad(1)=seprad/(1.d0+10.d0**(delmag(2,1)*0.4d0)) bseprad(2)=seprad-bseprad(1) bparad(1)=parad+pi bparad(2)=parad C DO IG=1,NFILT C VISTOT = CMPLX(0.d0,0.d0) inlambda=lambda(id,ig) color=ccoeff(3)*inlambda*inlambda+ & ccoeff(2)*inlambda + & ccoeff(1) unm=u(id,ig)*inlambda vnm=v(id,ig)*inlambda uvdist=dsqrt(unm**2+vnm**2) idf=int(dlambda(id,ig)) ! filter id DO IC = 1, NCOMP diamrad=(dcoeff(ic,3)*inlambda*inlambda+ & dcoeff(ic,2)*inlambda + & dcoeff(ic,1))*mas2rad C----------------------------------------------------------------------- C The following code is former subroutine BWVIS C C Calculate ul at wavelength lambda ul=lcoeff(ic,3)*inlambda*inlambda+ & lcoeff(ic,2)*inlambda + & lcoeff(ic,1) alpha=1.d0-ul beta=ul VIS = CMPLX(0.0D0, 0.0D0) C Do integration over bandpass i=1 cfactor=0.d0 do while(ftcurve(idf,i,2).ne.0) flambda=ftcurve(idf,i,1) CPHASE = TWOPI * BSEPRAD(IC) & * (UNM*DSIN(BPARAD(IC))+VNM*DCOS(BPARAD(IC))) & / flambda ARG = PI * DIAMRAD * UVDIST / flambda if(abs(arg).gt.1.d-14)then C DVISMOD = 2.0D0 * BESSJ1(ARG)/ARG ! Uniform disc C DVISMOD=EXP(-0.3D0*ARG**2) ! Gaussian disc dvismod = (alpha*bessj1(arg)/arg+beta*dsqrt(pi/2.D0)* & dsqrt(2.D0/(pi*arg))*(dsin(arg)/arg-dcos(arg))/ & dsqrt(arg*arg*arg))/(alpha/2.D0+beta/3.D0) else dvismod=1.d0 endif DVIS = CMPLX(DVISMOD*DCOS(CPHASE), DVISMOD*DSIN(CPHASE)) pcurve=(3.742d15/flambda**5)/ & (exp(1.439d0/(flambda*1.d-7*temp(ic)))-1.d0) VIS = VIS + DVIS * ftcurve(idf,i,2)/100.d0*pcurve CFACTOR = CFACTOR + ftcurve(idf,i,2)/100.d0*pcurve i=i+1 enddo VIS = VIS / CFACTOR C----------------------------------------------------------------------- c ...... Adjust binary component fluxes if(ic.ne.1)then flux=1.D0/(1.D0+10.D0**(-0.4D0*(delmag(ic,1)+color))) & *(10.D0**(-0.4D0*(delmag(ic,1)+color))) else flux=1.D0/(1.D0+10.D0**(-0.4D0*(delmag(2,1)+color))) endif VISTOT = VISTOT + FLUX*VIS END DO modvis2(id,ig) = DBLE(VISTOT)**2 + DIMAG(VISTOT)**2 modphas(id,ig) = DATAN2( DIMAG(VISTOT),DBLE(VISTOT) ) ENDDO c Restore input values for seprad and parad mseprad=seprado mparad=parado ENDDO c first=.false. RETURN c----------------------------------------------------------------------- c This code is to connect Bob Wilson's IMLC program to BINAFIT c All the parameters are set up for the case of Capella! 10 continue DO ID=1,NDATA DO IG=1,NFILT IDAY = MOD(DATE(ID), 100) IMON = MOD(DATE(ID)/100, 100) IYEAR = MOD(DATE(ID)/10000, 100) + 1900 DJUL = JD ( IDAY, IMON, IYEAR ) - 2440000.D0 DJUL = DJUL + HOURS(ID)/24.D0 c c IMLC parameters ifrad=0 ! controls output format (unused) nref=1 ! number of reflections if mref = 2 mref=1 ! simple treatment of reflection effect ifsmv1=1 ! star spot control #1 (unused) ifsmv2=0 ! star spot control #2 (unused) nsp1=0 nsp2=0 c if(first)then c idum=-1 c x=ran1(idum) c do i=1,nsp1 c spotparms(1,1,i)=90.+(ran1(idum)-0.5)*180. c spotparms(2,1,i)=ran1(idum)*360. c spotparms(3,1,i)=20. c spotparms(4,1,i)=0.8 c enddo c endif c do i=1,1 c spotparms(1,1,i)=110. c spotparms(2,1,i)=real(i-1)*45.+140. c spotparms(3,1,i)=40. c spotparms(4,1,i)=0.9 c enddo c spotparms(1,1,1)=180. spotparms(2,1,1)=0. spotparms(3,1,1)=20. spotparms(4,1,1)=0.8 c icor1=1 ! unknown icor2=1 ! unknown ld=1 ! linear cosine law for limb darkening modeb=2 ! detached binary ipb=0 ! no effect if mode = 0 ifat1=0 ! 0 = Blackbody radiation, 1 = stellar atmosphere ifat2=0 ! see ifat1 n1=10 ! number of latitude rows per hemisphere n2=10 ! see n1; longitude elements scales with sine of latitude bperiod=real(period) ! orbital period in days the=0. ! semi-duration of primary eclipse (unused) vunit=100. ! unit for radial velocity input and output (km/s) phstrt=0. ! start phase phstop=-1. ! stop phase (do only phase phn) phin=0.01 ! phase increment (we do only phn, however) e=real(e_true) ! eccentricity per=real(arg_peri) ! in degrees k1=27.40 ! km/s (hotter component of Capella) k2=26.05 ! km/s (cooler component) a1sini=13751.0*sqrt(1.d0-e**2)*k1*bperiod a2sini=13751.0*sqrt(1.d0-e**2)*k2*bperiod rsun=6.96e5 ! km (radius of sun) asemima=(a1sini+a2sini)/real(sin(inclin*degrad))/rsun f1=11.8 ! ratio of axial rotation rate to orbital rate f2= 1.0 vga=0.292 ! systemic velocity (in units of vunit) pshift=0. ! constant phase shift xincl=real(inclin) ! orbital inclination in degrees gr1=0.3 ! exponents in gravity darkening law gr2=0.3 ! see gr1, values for convective envelopes tavh=0.5700 ! effective temperature in 10000 K tavc=0.5100 ! see tavh alb1=0.5 ! bolometric albedo alb2=0.5 ! see alb1 ic=1 radius=real(dcoeff(ic,3)*lambda(id,ig)*lambda(id,ig)+ & dcoeff(ic,2)*lambda(id,ig) + & dcoeff(ic,1))/2. poth=real(a_true)/radius ! omega pot.= sem maj axis/rad of star ic=2 radius=real(dcoeff(ic,3)*lambda(id,ig)*lambda(id,ig)+ & dcoeff(ic,2)*lambda(id,ig) + & dcoeff(ic,1))/2. potc=real(a_true)/radius ! omega potential = semi major axis rm=k1/k2 ! mass ratio xbol1=0.4 ! coefficient in bolometric limb darkening law xbol2=0.4 ! see xbol1 ybol1=0.0 ! see xbol1, coefficient for logarithmic term ybol2=0.0 ! see ybol1, (unused if ld=1) wl=real(lambda(id,ig))/1000. ! wavlength in microns hlum=10.0 ! monochromatic luminosity, we will adjust this later clum=10.0 ! see hlum ic=1 ! monochromatic linear limb darkening coefficients xh=real(lcoeff(ic,3)*lambda(id,ig)*lambda(id,ig)+ & lcoeff(ic,2)*lambda(id,ig) + & lcoeff(ic,1)) ic=2 xc=real(lcoeff(ic,3)*lambda(id,ig)*lambda(id,ig)+ & lcoeff(ic,2)*lambda(id,ig) + & lcoeff(ic,1)) yh=0. ! see xh, coefficient for logarithmic term yc=0. ! see yh, (unused if ld=1) el3=0. ! third light (unused) zero=0. ! zero point for output magnitudes factor=1. ! scaling factor for normalized light column c c In IMLC, phase = 0 when eclipse occurs and arg_periastron = 90. c Therefore, we borrowed the following code from MODLOG in FNL c to calculate phase of the periastron, which is our zero point. c DTR=1.745329E-2 PERR=DTR*PER IF(E.NE.0.) GOTO 60 PERR=1.570796 60 CONTINUE TRC=1.570796-PERR 39 if(TRC.LT.0.) TRC=TRC+6.283185 if(trc.lt.0.) goto 39 40 if(trc.ge.6.283185) trc=trc-6.283185 if(trc.ge.6.283185) goto 40 HTRC=.5*TRC IF(ABS(1.570796-HTRC).LT.7.E-6) GOTO 101 IF(ABS(4.712389-HTRC).LT.7.E-6) GOTO 101 ECAN=2.*ATAN(SQRT((1.-E)/(1.+E))*TAN(HTRC)) GOTO 103 101 ECAN=3.141593 103 XMC=ECAN-E*SIN(ECAN) IF(XMC.LT.0.) XMC=XMC+6.283185 PHPER=1.-XMC/6.283185 PCONJ=(XMC+PERR)/6.283185-.25+PSHIFT 38 if(pconj.ge.1.) pconj=pconj-1. if(pconj.ge.1.) goto 38 41 if(pconj.lt.0.) pconj=pconj+1. if(pconj.lt.0.) goto 41 PHPERI=PHPER+PCONJ ! phase of periastron c Calculate normalized phase c phn=real(mod(djul-epoch,period)/period) ! phase of binary c if(e.eq.0.)then c phn=mod(phn+(per-90.)/360.,1.) c else c phn=mod(phn+phperi,1.) c endif c if(phn.lt.0.0)phn=phn+1. c if(phn.gt.1.0)phn=phn-1. c c Alternatively, calculate phase as is phn=real((djul-epoch)/period) ! phase of binary if(e.eq.0.)then phn=phn+(per-90.)/360. else phn=phn+phperi endif c c Before we call IMLC, we have to copy the parameters c into array lcparms c lcparms(1) = ifrad lcparms(2) = nref lcparms(3) = mref lcparms(4) = ifsmv1 lcparms(5) = ifsmv2 lcparms(6) = icor1 lcparms(7) = icor2 lcparms(8) = ld lcparms(9) = MODEB lcparms(10)= IPB lcparms(11)= IFAT1 lcparms(12)= IFAT2 lcparms(13)= N1 lcparms(14)= N2 lcparms(15)= BPERIOD lcparms(16)= THE lcparms(17)= VUNIT lcparms(18)= PHN lcparms(19)= PHSTRT lcparms(20)= PHSTOP lcparms(21)= PHIN lcparms(22)= E lcparms(23)= PER lcparms(24)= ASEMIMA lcparms(25)= F1 lcparms(26)= F2 lcparms(27)= VGA lcparms(28)= PSHIFT lcparms(29)= XINCL lcparms(30)= GR1 lcparms(31)= GR2 lcparms(32)= tavh lcparms(33)= tavc lcparms(34)= alb1 lcparms(35)= alb2 lcparms(36)= poth lcparms(37)= potc lcparms(38)= rm lcparms(39)= xbol1 lcparms(40)= xbol2 lcparms(41)= ybol1 lcparms(42)= ybol2 lcparms(43)= WL lcparms(44)= HLUM lcparms(45)= CLUM lcparms(46)= XH lcparms(47)= xc lcparms(48)= yh lcparms(49)= yc lcparms(50)= EL3 lcparms(51)= ZERO lcparms(52)= FACTOR lcparms(53)= NSP1 lcparms(54)= NSP2 do l=1,54 if(lcparms(l).ne.lcoparms(l))then call imlc(lcparms,spotparms) ! call only if parameters changed goto 20 endif enddo goto 30 20 if(ilc.gt.maxlc)then write(outc,*)'***ERROR: maximum number of grid points exceeded!' stop endif do l=1,54 lcoparms(l)=lcparms(l) enddo c c Determine total flux per component, c scale coordinates and rotate counterclockwise by asc_node-90 30 rangle=(asc_node-90.d0)*degrad do ic=1,2 flctot(ic)=0.d0 enddo do i=1,ilc flctot(clc(i))=flctot(clc(i))+flc(i) hxlc(i)=(-ylc(i)*sin(rangle)+xlc(i)*cos(rangle))*a_true*mas2rad hylc(i)=(+xlc(i)*sin(rangle)+ylc(i)*cos(rangle))*a_true*mas2rad enddo color=ccoeff(3)*lambda(id,ig)*lambda(id,ig)+ & ccoeff(2)*lambda(id,ig) + & ccoeff(1) c 1st component flux=1.D0/(1.D0+10.D0**(-0.4D0*(delmag(2,1)+color))) flctot(1)=flux/flctot(1) c 2nd component flux=1.D0/(1.D0+10.D0**(-0.4D0*(delmag(2,1)+color))) & *(10.D0**(-0.4D0*(delmag(2,1)+color))) flctot(2)=flux/flctot(2) c c Do Fourier transformation, including bandwidth smearing c vistot=cmplx(0.d0,0.d0) idf=int(dlambda(id,ig)) ! find filter id unm=u(id,ig)*lambda(id,ig) vnm=v(id,ig)*lambda(id,ig) l=1 nfactor=0.d0 do while(ftcurve(idf,l,2).ne.0) ulam=unm/ftcurve(idf,l,1) vlam=vnm/ftcurve(idf,l,1) vis=cmplx(0.d0,0.d0) do i=1,ilc vis=vis+flctot(clc(i))*flc(i) & *cmplx(cos(2.d0*pi*(-ulam*hxlc(i)+vlam*hylc(i))), & sin(2.d0*pi*(-ulam*hxlc(i)+vlam*hylc(i)))) enddo vistot=vistot+vis*ftcurve(idf,l,2)/100.d0 nfactor=nfactor +ftcurve(idf,l,2)/100.d0 vistot=vistot+vis nfactor=nfactor+1.d0 l=l+1 enddo VISTOT = VISTOT / NFACTOR modvis2(id,ig) = DBLE(VISTOT)**2 + DIMAG(VISTOT)**2 modphas(id,ig) = DATAN2( DIMAG(VISTOT),DBLE(VISTOT) ) ENDDO ENDDO first=.false. return END