SUBROUTINE BWVIS(IC,ID,IG,VISIB) C This program calculates the visibility (NOT visibility squared!) C of a uniform disk source not at the phase center, observed C with non-negligible bandwidth. C C Version 1.0 12 Dec. 1990 J.T. Armstrong C C The integral can be written as: C V = N * Int { dnu C Int [ sigma dsigma C exp(-2*pi*i*sep*uvdist*(nu/c)*cos(posang - uvangle)) C Int (dgamma C exp(-2*pi*i*sigma*uvdist1*(nu/c)*cos(gamma)) ) ] } C where: C N is the normalization (see below) C nu is frequency (Hz) C c is the speed of light (m/s) C The integral is expressed in two systems of cylindrical coordinates. C The first is centered at the phase center. The quantities expressed C in this system are: C uvdist (m): the length of the baseline component in the plane of C the sky, i.e., sqrt(u^2 + v^2), but in meters C uvangle(rad): the position angle of the projected baseline, i.e., C atan(u/v) C sep (mas): the distance of the source from the phase center C posang (rad): the position angle of the source C The second cylindrical coordinate system is centered on the source. C The quantities expressed in this system are: C uvdist1 (m): the length of the baseline component in the plane of C the sky, i.e., sqrt(u^2 + v^2), but in meters, and C measured at the source center. This will be very C little different from uvdist. C sigma (mas): the distance of a source point from the center of C the source C gamma (rad): the position angle of a source point C Other quantities of interest are: C lambda0 (nm): the center wavelength C dlambda (nm): the filter width C D (mas): the diameter of the source C N: the normalization, equal to ((nu(hi)-nu(lo))*pi*r^2)^-1. C C Notice that the gamma integral results in C 2*pi*BESJ0(2*pi*sigma*uvdist1*(nu/c)), where BESJ0() is the Bessel C function of the first kind and zeroth order. C C NOTE: I assume a square bandpass centered in wavelength on lambda0. C This leads to various complications in the frequency integral limits. C I also assume a uniform brightness across the source. Either or both C of these could be changed by introducing a function T(nu) for the C filter transmission and/or I(sigma/R) for the brightness as a function C of radius. The normalization would have to be adjusted accordingly; C in particular, C N = 2*pi*Int { dnu T(nu) C Int [sigma dsigma I(sigma/R) C Int (dgamma) ] }, C where the 2*pi comes from integrating over gamma. C More complicated source shapes, such as an elliptical star, would C require doing the sigma integral inside the gamma integral. INCLUDE 'BINFIT.INC' C ARRAYS: REAL*8 SIGMA(SIGMAX), GAMMA(GAMMAX) C U-V PLANE QUANTITIES: REAL*8 UVDIST, UVDIST1, UVANGLE C FREQUENCY: REAL*8 NU C PHASES: REAL*8 CPHASE, BESJ0ARG, BESJ1ARG, VISPHA C INTEGRANDS AND VISIBILITIES: REAL*8 BESSJ1, DVISMOD, VISMOD, INTMOD COMPLEX*16 INTEGRAND, VISIB, DVIS, DINT C INDEXES: INTEGER*4 ISIGMA, INU, IGAMMA C Used for limb-darkening code logical first real*8 ul(mxmod,mxfilt) common/limb/ul real*8 arg,alpha,beta integer*4 iic,iig c save seprad, parad, since their new values due to orbital motion c is used real*8 seprado,parado data first/.true./ 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 c if(first)then first=.false. write(outc,*)'Please enter limb darkening coefficients!' do iic=1,ncomp do iig=1,nfilt ul(iic,iig)=0.D0 write(outc,'(a,i1,a,i1,a,$)')' Component ',iic, & ' Filter ',iig,': [0] ' read(inc,*)ul(iic,iig) enddo enddo endif C Set up arrays of integration variables: C==================================== C To be used for limb-darkened stars: C DO ISIGMA = 1, SIGSTEPS C SIGMA(ISIGMA) = (ISIGMA-1)*(D/2.0)/(SIGSTEPS-1) C END DO C INTEGRAND = CMPLX(0.0D0,0.0D0) C DO ISIGMA = 1, SIGSTEPS C ARG = TWOPI * SIGMA(ISIGMA) * UVDIST1 * NU(INU)/C C INTMOD = TWOPI * BESJ0(BESJ0ARG) C INTEGRAND = INTEGRAND + INTMOD C END DO C==================================== C==================================== C To be used for noncircular sources: C DO IGAMMA = 1, GAMSTEPS C GAMMA(IGAMMA) = TWOPI*(IGAMMA-1)/GAMSTEPS C END DO C==================================== C Calculate the UV distance and position angle. D WRITE(OUTC,'(A,3(1X,I4))') D 1 ' BWVIS: comp., scan, filter: ', D 2 IC,ID,IG C WRITE(OUTC,'(A,2(1X,G12.5))') ' U, V for this scan, filter: ', C 1 UDATA(ID,IG),VDATA(ID,IG) UVDIST = DSQRT(UDATA(ID,IG)**2 + VDATA(ID,IG)**2) UVDIST1 = UVDIST ! Diff. should be very small UVANGLE = DATAN2(VDATA(ID,IG),UDATA(ID,IG)) D WRITE(OUTC,'(A,G10.5,2X,G10.5,2X)') D 1 ' UVDIST, UVANGLE: ', UVDIST, UVANGLE alpha=1.D0-ul(ic,ig) beta=ul(ic,ig) C Do the integration: c write(15,*)'alpha=',alpha,' beta=',beta,' lambda=',c/nu0(ig)*1.d9 c write(15,*)'u=',udata(id,ig),' v=',vdata(id,ig) c write(15,*)'r=',seprad(ic,ig)/mas2rad,'theta=',parad(ic,ig)/degrad c write(15,*)'diam=',diamrad(ic,ig)/mas2rad VISIB = CMPLX(0.0D0, 0.0D0) NU = NULO(IG) - 0.5 * DNU(IG) D WRITE(OUTC,'(A,G12.5)') ' NU = ',NU D WRITE(OUTC,*) ' Start frequency integration' DO INU = 1, NUSTEPS NU = NU + DNU(IG) CPHASE = TWOPI * SEPRAD(IC,IG) * (NU/NU0(IG)) 1 * ( UDATA(ID,IG) * DSIN(PARAD(IC,IG)) 1 + VDATA(ID,IG) * DCOS(PARAD(IC,IG)) ) BESJ1ARG = PI * DIAMRAD(IC,IG) * UVDIST1 * (NU/NU0(IG)) D WRITE(OUTC,'(A,G14.6)') ' Arg. of J1: ',BESJ1ARG 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), DVISMOD*DSIN(CPHASE)) VISIB = VISIB + DVIS END DO D WRITE(OUTC,*) ' Finish frequency integration' VISIB = VISIB / NUSTEPS D WRITE(OUTC,'(A,I1,A,I1,A,5(G11.4))') D 1 ' D,ANG,S,PA,VISIB (1,', D 1 IC,',',IG,': ',UVDIST,UVANGLE,SEPRAD(IC,IG),PARAD(IC,IG),VISIB seprad(ic,ig)=seprado parad(ic,ig)=parado RETURN END