vis_calc_ext.mws

VLTI Baseline Vector and OPD

(16.08.2007)

1. Initialisations

> restart:

> Digits := 16;

Digits := 16

> deg2rad := evalf(180./Pi);

deg2rad := 57.29577951308233

> mas2rad := deg2rad*1000.*3600.;

mas2rad := 206264806.2470964

> with(linalg): # Linear algebra package

Warning, the protected names norm and trace have been redefined and unprotected

> with(inttrans): # Integral transforms

Warning, the name hilbert has been redefined

> latitude := -24.62528; # Latitude of Paranal obs., in degrees (about 24 degrees South, South counts as negative)

latitude := -24.62528

> longitude := -70.40278; # Longitude of Paranal obs., in degrees (about 70 degrees West, West counts as negative)

longitude := -70.40278

2. Computing the number of days since J2000.0 (Jan 1st 2000, 12:00 UT)

> # For the definition of ydays we voluntarily avoid the leap year 2004

> ydays := yy -> piecewise(yy=1999,-366.5,yy=2000,-1.5,yy=2001,364.5,yy=2002,729.5,yy=2003,1094.5);

ydays := proc (yy) options operator, arrow; piecewi...

> # On leap years add one to the days of the month starting with March

> mdays := mm -> piecewise(mm=1,0,mm=2,31,mm=3,59,mm=4,90,mm=5,120,mm=6,151,mm=7,181,mm=8,212,mm=9,243,mm=10,273,mm=11,304,mm=12,334);

mdays := proc (mm) options operator, arrow; piecewi...

> yy := 2003: # 2003-01-01

> mm := 01:

> dd := 01:

>

> # In this application we leave UT undetermined for the moment

> # UT := 23.00 + 20/60. + 44.5/3600.; # Universal Time of observation (23h20min44.500 for alfpsa)

>

> days2K := ydays(yy) + mdays(mm) + dd + UT/24.0; # Number of days since J2000.0 (J2000.0 corresponds to 1200 hrs UT on Jan 1st 2000 AD)

days2K := 1095.5+.4166666666666667e-1*UT

3. Computing the Sun coordinates (right ascension, declination)

> # Computes Sun coordinates

> Ls := 280.461 + 0.9856474 * days2K; # Mean longitude of the sun

Ls := 1360.23772670+.4106864166666667e-1*UT

> Ls := Ls - 360.*floor(Ls/360.); # Ls (in degrees) is in the interval [0.,360.]

Ls := 1360.23772670+.4106864166666667e-1*UT-360.*fl...

> gs := 357.528 + 0.9856003 * days2K; # Mean anomaly of the sun

gs := 1437.25312865+.4106667916666667e-1*UT

> gs := gs - 360.*floor(gs/360.); # gs (in degrees) is in the interval [0.,360.]

gs := 1437.25312865+.4106667916666667e-1*UT-360.*fl...

> lambda0 := Ls + 1.915 * sin(gs/deg2rad) + 0.020 * sin(2*gs/deg2rad); # Ecliptic longitude of the Sun, in degrees

lambda0 := 1360.23772670+.4106864166666667e-1*UT-36...
lambda0 := 1360.23772670+.4106864166666667e-1*UT-36...
lambda0 := 1360.23772670+.4106864166666667e-1*UT-36...

> lambda0 := lambda0 - 360.*floor(lambda0/360.);

lambda0 := 1360.23772670+.4106864166666667e-1*UT+1....
lambda0 := 1360.23772670+.4106864166666667e-1*UT+1....
lambda0 := 1360.23772670+.4106864166666667e-1*UT+1....
lambda0 := 1360.23772670+.4106864166666667e-1*UT+1....

> epsilon0 := 23.439 - 0.0000004 * days2K; # Obliquity of the ecliptic plane

epsilon0 := 23.43856180-.1666666666666667e-7*UT

> Ys := cos(epsilon0/deg2rad) * sin(lambda0/deg2rad);

Ys := cos(-.4090800753421685+.2908882086657216e-9*U...
Ys := cos(-.4090800753421685+.2908882086657216e-9*U...
Ys := cos(-.4090800753421685+.2908882086657216e-9*U...
Ys := cos(-.4090800753421685+.2908882086657216e-9*U...
Ys := cos(-.4090800753421685+.2908882086657216e-9*U...

> Xs := cos(lambda0/deg2rad);

Xs := cos(23.74062694075778+.7167830164050647e-3*UT...
Xs := cos(23.74062694075778+.7167830164050647e-3*UT...
Xs := cos(23.74062694075778+.7167830164050647e-3*UT...
Xs := cos(23.74062694075778+.7167830164050647e-3*UT...

> alpha0 := evalf(arctan(Ys,Xs)*deg2rad);

alpha0 := 57.29577951308233*arctan(cos(-.4090800753...
alpha0 := 57.29577951308233*arctan(cos(-.4090800753...
alpha0 := 57.29577951308233*arctan(cos(-.4090800753...
alpha0 := 57.29577951308233*arctan(cos(-.4090800753...
alpha0 := 57.29577951308233*arctan(cos(-.4090800753...
alpha0 := 57.29577951308233*arctan(cos(-.4090800753...
alpha0 := 57.29577951308233*arctan(cos(-.4090800753...
alpha0 := 57.29577951308233*arctan(cos(-.4090800753...
alpha0 := 57.29577951308233*arctan(cos(-.4090800753...

> delta0 := deg2rad*arcsin(sin(epsilon0/deg2rad)*sin(lambda0/deg2rad));

delta0 := -57.29577951308233*arcsin(sin(-.409080075...
delta0 := -57.29577951308233*arcsin(sin(-.409080075...
delta0 := -57.29577951308233*arcsin(sin(-.409080075...
delta0 := -57.29577951308233*arcsin(sin(-.409080075...
delta0 := -57.29577951308233*arcsin(sin(-.409080075...

4. Computing the Local Sidereal Time (LST)

> # Compute Local Sidereal Time (LST)

> LST := 100.46 + 0.985647 * days2K + longitude + 15.0*UT; # LST is expressed in degrees

LST := 1109.8335085+15.04106862500000*UT

> LST := LST-3*360.;

LST := 29.8335085+15.04106862500000*UT

> LSTS := LST*240.; # LST in seconds, as given in the FITS header

LSTS := 7160.0420400+3609.856470000000*UT

>

> LSTH := (LST)/15;

LSTH := 1.988900566666667+1.002737908333333*UT

> LSTM := (LSTH-20.)*60;

LSTM := -1080.665966000000+60.16427449999998*UT

> LSTS := (LSTM-36.)*60.; # LST is 00:40:08.0745

LSTS := -66999.95796000000+3609.856469999999*UT

5. Sun altitude and azimuth

> # Coumpute Hour Angle of the Sun, Altitude and Azimuth (HAS, ALTS, AZS)

> HAS := LST - alpha0;

HAS := 29.8335085+15.04106862500000*UT-57.295779513...
HAS := 29.8335085+15.04106862500000*UT-57.295779513...
HAS := 29.8335085+15.04106862500000*UT-57.295779513...
HAS := 29.8335085+15.04106862500000*UT-57.295779513...
HAS := 29.8335085+15.04106862500000*UT-57.295779513...
HAS := 29.8335085+15.04106862500000*UT-57.295779513...
HAS := 29.8335085+15.04106862500000*UT-57.295779513...
HAS := 29.8335085+15.04106862500000*UT-57.295779513...
HAS := 29.8335085+15.04106862500000*UT-57.295779513...

> sas := sin(delta0/deg2rad)*sin(latitude/deg2rad)+cos(delta0/deg2rad)*cos(latitude/deg2rad)*cos(HAS/deg2rad);

sas := .4166819241881570*sin(1.000000000000000*arcs...
sas := .4166819241881570*sin(1.000000000000000*arcs...
sas := .4166819241881570*sin(1.000000000000000*arcs...
sas := .4166819241881570*sin(1.000000000000000*arcs...
sas := .4166819241881570*sin(1.000000000000000*arcs...
sas := .4166819241881570*sin(1.000000000000000*arcs...
sas := .4166819241881570*sin(1.000000000000000*arcs...
sas := .4166819241881570*sin(1.000000000000000*arcs...
sas := .4166819241881570*sin(1.000000000000000*arcs...
sas := .4166819241881570*sin(1.000000000000000*arcs...
sas := .4166819241881570*sin(1.000000000000000*arcs...
sas := .4166819241881570*sin(1.000000000000000*arcs...
sas := .4166819241881570*sin(1.000000000000000*arcs...
sas := .4166819241881570*sin(1.000000000000000*arcs...
sas := .4166819241881570*sin(1.000000000000000*arcs...
sas := .4166819241881570*sin(1.000000000000000*arcs...
sas := .4166819241881570*sin(1.000000000000000*arcs...
sas := .4166819241881570*sin(1.000000000000000*arcs...
sas := .4166819241881570*sin(1.000000000000000*arcs...

> ALTS := arcsin(sas);

ALTS := arcsin(.4166819241881570*sin(1.000000000000...
ALTS := arcsin(.4166819241881570*sin(1.000000000000...
ALTS := arcsin(.4166819241881570*sin(1.000000000000...
ALTS := arcsin(.4166819241881570*sin(1.000000000000...
ALTS := arcsin(.4166819241881570*sin(1.000000000000...
ALTS := arcsin(.4166819241881570*sin(1.000000000000...
ALTS := arcsin(.4166819241881570*sin(1.000000000000...
ALTS := arcsin(.4166819241881570*sin(1.000000000000...
ALTS := arcsin(.4166819241881570*sin(1.000000000000...
ALTS := arcsin(.4166819241881570*sin(1.000000000000...
ALTS := arcsin(.4166819241881570*sin(1.000000000000...
ALTS := arcsin(.4166819241881570*sin(1.000000000000...
ALTS := arcsin(.4166819241881570*sin(1.000000000000...
ALTS := arcsin(.4166819241881570*sin(1.000000000000...
ALTS := arcsin(.4166819241881570*sin(1.000000000000...
ALTS := arcsin(.4166819241881570*sin(1.000000000000...
ALTS := arcsin(.4166819241881570*sin(1.000000000000...
ALTS := arcsin(.4166819241881570*sin(1.000000000000...
ALTS := arcsin(.4166819241881570*sin(1.000000000000...

> altsdeg := ALTS*deg2rad;

altsdeg := 57.29577951308233*arcsin(.41668192418815...
altsdeg := 57.29577951308233*arcsin(.41668192418815...
altsdeg := 57.29577951308233*arcsin(.41668192418815...
altsdeg := 57.29577951308233*arcsin(.41668192418815...
altsdeg := 57.29577951308233*arcsin(.41668192418815...
altsdeg := 57.29577951308233*arcsin(.41668192418815...
altsdeg := 57.29577951308233*arcsin(.41668192418815...
altsdeg := 57.29577951308233*arcsin(.41668192418815...
altsdeg := 57.29577951308233*arcsin(.41668192418815...
altsdeg := 57.29577951308233*arcsin(.41668192418815...
altsdeg := 57.29577951308233*arcsin(.41668192418815...
altsdeg := 57.29577951308233*arcsin(.41668192418815...
altsdeg := 57.29577951308233*arcsin(.41668192418815...
altsdeg := 57.29577951308233*arcsin(.41668192418815...
altsdeg := 57.29577951308233*arcsin(.41668192418815...
altsdeg := 57.29577951308233*arcsin(.41668192418815...
altsdeg := 57.29577951308233*arcsin(.41668192418815...
altsdeg := 57.29577951308233*arcsin(.41668192418815...
altsdeg := 57.29577951308233*arcsin(.41668192418815...

> plot(altsdeg,UT=0..24);

[Maple Plot]

> # Let's consider night time is from 00:00 UT to 09:00 UT

6. Target altitude and azimuth

> # Compute Hour Angle of the target

> # Target is Sirius

> rah := 06.00: # RA = 6 hours and...

> ram := 45.00: # 45 minutes and ...

> ras := 08.90: # 8.9 seconds

> alpha := rah*15+ram*15/60+ras*15/3600; # Target right ascension, in degrees

alpha := 101.2870833333333

>

> decd := -16: # Declination is -16 degrees and ..

> decm := 42: # 42 minutes and...

> decs := 58.5: # 58.5 seconds

> dd := abs(decd); # we keep the absolute value and sign of the declination

dd := 16

> ds := sign(decd);

ds := -1

> delta := ds*(dd+decm/60.+decs/3600.); # Target declination, in degrees

delta := -16.71625000000000

> HA := LST - alpha; # Hour Angle in degrees

HA := -71.4535748333333+15.04106862500000*UT

> HA := HA + 360.;

HA := 288.5464251666667+15.04106862500000*UT

> HA/15; # Hour Angle in Hours

19.23642834444445+1.002737908333333*UT

> plot(HA/15,UT=5..9); # Plot hour angle in hours: 24 hours=360 degrees

[Maple Plot]

> # Compute Altitude and Azimuth (ALT, AZ) of the target

> sa := sin(delta/deg2rad)*sin(latitude/deg2rad)+cos(delta/deg2rad)*cos(latitude/deg2rad)*cos(HA/deg2rad);

sa := .1198511227525249+.8706366664797497*cos(5.036...

> altitude := arcsin(sa);

altitude := arcsin(.1198511227525249+.8706366664797...

> altdeg := altitude*deg2rad; # Force floating point evaluation of altitude, convert to degrees

altdeg := 57.29577951308233*arcsin(.119851122752524...

> plot(altdeg,UT=0..24);

[Maple Plot]

> # During the night time (00:00 UT to 09:00 UT), the target is at reasonable elevation (above 20 degrees) only until 07:00 UT. We will later on restrict our time exploration to the range 00:00 - 07:00 UT.

> tazx := cos(delta/deg2rad)*sin(HA/deg2rad);

tazx := .9577409562815381*sin(5.036085164017760+.26...

> tazy := sin(latitude/deg2rad)*cos(delta/deg2rad)*cos(HA/deg2rad)-cos(latitude/deg2rad)*sin(delta/deg2rad);

tazy := -.3990733445371968*cos(5.036085164017760+.2...

> azimuth := -arctan(tazy,tazx); # -atan2(tazy,tazx);

azimuth := -arctan(-.3990733445371968*cos(5.0360851...

> evalf(arctan(-1,1)*deg2rad);

-45.00000000000000

> azd := azimuth*deg2rad;

azd := -57.29577951308233*arctan(-.3990733445371968...

> azdeg := piecewise(azd<0,azd+360,azd>=0,azd);

azdeg := PIECEWISE([-57.29577951308233*arctan(-.399...
azdeg := PIECEWISE([-57.29577951308233*arctan(-.399...
azdeg := PIECEWISE([-57.29577951308233*arctan(-.399...
azdeg := PIECEWISE([-57.29577951308233*arctan(-.399...

> plot(azdeg,UT=0..7);

[Maple Plot]

> plot(270-azdeg,UT=0..10); # VisCalc azimuth plot

[Maple Plot]

8. Compute baseline vector.

> nu := -18.984:

> psi := evalf(nu/deg2rad); # Rotation of the platform (can be set to nu=-18.984 degrees as above)

psi := -.3313333051986035

> R := matrix(3,3,[[cos(psi),sin(psi),0],[-sin(psi),cos(psi),0],[0,0,1]]);

R := matrix([[.9456094545111512, -.3253041031698232...

> # Exact coordinates for the stations are U1 -16.000 -16.000 -9.925 -20.335, U3 64.0013 47.9725 44.915 66.183

> ntel1 := vector([-16,-16,0]); # Nominal U,V position for station UT1

ntel1 := vector([-16, -16, 0])

> ntel2 := vector([64.0013,47.9725,0]); # Nominal U,V position for station UT3

ntel2 := vector([64.0013, 47.9725, 0])

> rtel1 := evalm(R &* ntel1); # N-E coordinates of station E0, as given in the FITS header

rtel1 := vector([-9.924885621461249, -20.3346169228...

> rtel2 := evalm(R &* ntel2); # N-E coordinates of station G0, as given in the FITS header

rtel2 := vector([44.91458329169020, 66.183135054739...

>

>

> BN := rtel2[2] - rtel1[2]; # N-S ground projection of the baseline vector

BN := 86.51775197763460

> BE := rtel2[1] - rtel1[1]; # E-W ground projection of the baseline vector

BE := 54.83946891315145

> BA := rtel2[3] - rtel1[3]; # Elevation difference between telescope 1 and 2

BA := 0.

> Bg:= vector([BN,BE,BA]);

Bg := vector([86.51775197763460, 54.83946891315145,...

> harad := HA/deg2rad;

harad := 5.036085164017760+.2625161705246662*UT

> decrad := delta/deg2rad;

decrad := -.2917536010865021

> decrad := -.2917536010865021;

decrad := -.2917536010865021

> #harad := -15./60./deg2rad; # Overwrite, H = 1 min

> #decrad := 0./deg2rad; # Overwrite

> lat := latitude/deg2rad;

lat := -.4297922152255092

> L := matrix(3,3,[[-sin(lat),0,cos(lat)],[0,-1,0],[cos(lat),0,sin(lat)]]); # Belle, Eq. 4.108, p. 44

L := matrix([[.4166819241881570, 0, .90905234945785...

> B := evalm(L &* Bg);

B := vector([36.05038337047451, -54.83946891315145,...

> T := matrix(3,3,[[sin(harad),-cos(harad),0],[-sin(decrad)*cos(harad),-sin(decrad)*sin(harad),cos(decrad)],[cos(decrad)*cos(harad),cos(decrad)*sin(harad),sin(decrad)]]);

T := matrix([[sin(5.036085164017760+.26251617052466...

> UVW := evalm(T &* B);

UVW := vector([36.05038337047451*sin(5.036085164017...
UVW := vector([36.05038337047451*sin(5.036085164017...
UVW := vector([36.05038337047451*sin(5.036085164017...

> u := UVW[1];

u := 36.05038337047451*sin(5.036085164017760+.26251...

> plot(u,UT=0..9);

[Maple Plot]

> v:= UVW[2];

v := 10.36924971254399*cos(5.036085164017760+.26251...

> plot(v,UT=0..9);

[Maple Plot]

> plot([[u,v,UT=0..9],[-u,-v,UT=0..9]],-100..100,-100..100);

[Maple Plot]

>

> w := UVW[3];

w := 34.52692864355432*cos(5.036085164017760+.26251...

> plot(w,UT=3..3.5);

[Maple Plot]

> rho := sqrt(u*u+v*v);

rho := sqrt((36.05038337047451*sin(5.03608516401776...
rho := sqrt((36.05038337047451*sin(5.03608516401776...

> #plot(rho,decrad=-30/deg2rad..-20/deg2rad);

9. Model Uniform Disk

> lambda := 2.2e-6; # Central wavelength in meters (1 micron = 1e-6 meter)

lambda := .22e-5

> theta := 6.00/mas2rad; # 6 mas, converted to radian

theta := .2908882086657215e-7

> rho := sqrt(u*u + v*v);

rho := sqrt((36.05038337047451*sin(5.03608516401776...
rho := sqrt((36.05038337047451*sin(5.03608516401776...

> C(rho);

C(sqrt((36.05038337047451*sin(5.036085164017760+.26...
C(sqrt((36.05038337047451*sin(5.036085164017760+.26...

> plot(rho,UT=3..3.5);

[Maple Plot]

> uda := Pi*theta*abs(rho)/lambda;

uda := .1322219130298734e-1*Pi*sqrt(abs((36.0503833...
uda := .1322219130298734e-1*Pi*sqrt(abs((36.0503833...

> evalf(uda);

.4153873906182388e-1*sqrt(abs((36.05038337047451*si...
.4153873906182388e-1*sqrt(abs((36.05038337047451*si...

> udvs := 2*BesselJ(1,uda)/uda;

udvs := 151.2608579145374*BesselJ(1,.13222191302987...
udvs := 151.2608579145374*BesselJ(1,.13222191302987...
udvs := 151.2608579145374*BesselJ(1,.13222191302987...
udvs := 151.2608579145374*BesselJ(1,.13222191302987...

> udv := evalf(abs(udvs));

udv := 48.14782646683894*abs(BesselJ(1.,.4153873906...
udv := 48.14782646683894*abs(BesselJ(1.,.4153873906...
udv := 48.14782646683894*abs(BesselJ(1.,.4153873906...
udv := 48.14782646683894*abs(BesselJ(1.,.4153873906...

> plot(udv,UT=0..10); # Visibility values

[Maple Plot]

> # Phase

> # if (udvs < 0.0) phase = pi; else if (udvs >= 0.0) phase = 0.0;

9. Model Gaussian

> lambda := 2.20e-6; # Central wavelength in meters (1 micron = 1e-6 meter)

lambda := .220e-5

> uCoord := u/lambda;

uCoord := 16386537.89567023*sin(5.036085164017760+....

> vCoord := v/lambda;

vCoord := 4713295.323883631*cos(5.036085164017760+....

> rho := sqrt(uCoord*uCoord + vCoord*vCoord);

rho := sqrt((16386537.89567023*sin(5.03608516401776...
rho := sqrt((16386537.89567023*sin(5.03608516401776...

rho := sqrt((16391383.02851202*sin(5.03608516401776...
rho := sqrt((16391383.02851202*sin(5.03608516401776...

> plot(evalf(rho),UT=0..1); # Baseline plot

[Maple Plot]

>

> gaussFwhm := 6.00/mas2rad; # 6 mas, converted to radian

gaussFwhm := .2908882086657215e-7

> udv := exp( (-2 * Pi * Pi * gaussFwhm * gaussFwhm * rho * rho)/(8 * log(2)) );

udv := exp(-.2115398748518808e-15*Pi^2*((16386537.8...
udv := exp(-.2115398748518808e-15*Pi^2*((16386537.8...

> plot(evalf(udv),UT=0..10,0..0.025); # Visibility values

[Maple Plot]

> phase := 0; # for a Gaussian the phase is always zero

phase := 0

10. Model Binary Disc

> lambda := 2.2e-6; # Central wavelength in meters (1 micron = 1e-6 meter)

lambda := .22e-5

> uCoord := u/lambda;

uCoord := 16386537.89567023*sin(5.036085164017760+....

> vCoord := v/lambda;

vCoord := 4713295.323883631*cos(5.036085164017760+....

> rho := sqrt(uCoord*uCoord + vCoord*vCoord);

rho := sqrt((16386537.89567023*sin(5.03608516401776...
rho := sqrt((16386537.89567023*sin(5.03608516401776...

> brightRatio := 1.0; # Brightness ratio

brightRatio := 1.0

> r1 := 4.0/2/mas2rad; # Radius 1 (half-diameter), converted to radian

r1 := .9696273622190718e-8

> r2 := 3.0/2/mas2rad; # Radius 2 (half-diameter), converted to radian

r2 := .7272205216643038e-8

> theta := 37; # Position angle, converted to radian

theta := 37

> sep := 10/mas2rad; # Separation 12 mas, converted to radian

sep := .4848136811095359e-7

> ud1 := 2.0*Pi*rho*r1;

ud1 := .1939254724438144e-7*Pi*sqrt((16386537.89567...
ud1 := .1939254724438144e-7*Pi*sqrt((16386537.89567...

> ud2 := 2.0*Pi*rho*r2;

ud2 := .1454441043328608e-7*Pi*sqrt((16386537.89567...
ud2 := .1454441043328608e-7*Pi*sqrt((16386537.89567...

> bessel1 := BesselJ(1,ud1);

bessel1 := BesselJ(1,.1939254724438144e-7*Pi*sqrt((...
bessel1 := BesselJ(1,.1939254724438144e-7*Pi*sqrt((...

> bessel2 := BesselJ(1,ud2);

bessel2 := BesselJ(1,.1454441043328608e-7*Pi*sqrt((...
bessel2 := BesselJ(1,.1454441043328608e-7*Pi*sqrt((...

> udv_den := ( Pi*Pi* (brightRatio*brightRatio*r1*r1*r1*r1*uCoord*uCoord + brightRatio*brightRatio*r1*r1*r1*r1*vCoord*vCoord + 2*brightRatio*r1*r1*r2*r2*uCoord*uCoord + 2*brightRatio*r1*r1*r2*r2*vCoord*vCoord + r2*r2*r2*r2*uCoord*uCoord + r2*r2*r2*r2*vCoord*vCoord ));

udv_den := Pi^2*(.2158040058465925e-31*(16386537.89...
udv_den := Pi^2*(.2158040058465925e-31*(16386537.89...

> udv_num := 2*r2*bessel2*brightRatio*r1*bessel1*cos(2*Pi*sep*(cos(theta)*uCoord + sin(theta)*vCoord)) + brightRatio*brightRatio*r1*r1*bessel1*bessel1 + r2*r2*bessel2*bessel2;

udv_num := .1410265832345873e-15*BesselJ(1,.1454441...
udv_num := .1410265832345873e-15*BesselJ(1,.1454441...
udv_num := .1410265832345873e-15*BesselJ(1,.1454441...
udv_num := .1410265832345873e-15*BesselJ(1,.1454441...
udv_num := .1410265832345873e-15*BesselJ(1,.1454441...
udv_num := .1410265832345873e-15*BesselJ(1,.1454441...
udv_num := .1410265832345873e-15*BesselJ(1,.1454441...
udv_num := .1410265832345873e-15*BesselJ(1,.1454441...
udv_num := .1410265832345873e-15*BesselJ(1,.1454441...
udv_num := .1410265832345873e-15*BesselJ(1,.1454441...
udv_num := .1410265832345873e-15*BesselJ(1,.1454441...

>

> udv2 := udv_num / udv_den;

udv2 := (.1410265832345873e-15*BesselJ(1,.145444104...
udv2 := (.1410265832345873e-15*BesselJ(1,.145444104...
udv2 := (.1410265832345873e-15*BesselJ(1,.145444104...
udv2 := (.1410265832345873e-15*BesselJ(1,.145444104...
udv2 := (.1410265832345873e-15*BesselJ(1,.145444104...
udv2 := (.1410265832345873e-15*BesselJ(1,.145444104...
udv2 := (.1410265832345873e-15*BesselJ(1,.145444104...
udv2 := (.1410265832345873e-15*BesselJ(1,.145444104...
udv2 := (.1410265832345873e-15*BesselJ(1,.145444104...
udv2 := (.1410265832345873e-15*BesselJ(1,.145444104...
udv2 := (.1410265832345873e-15*BesselJ(1,.145444104...
udv2 := (.1410265832345873e-15*BesselJ(1,.145444104...
udv2 := (.1410265832345873e-15*BesselJ(1,.145444104...

> udv := sqrt(udv2);

udv := 1/Pi*((.1410265832345873e-15*BesselJ(1,.1454...
udv := 1/Pi*((.1410265832345873e-15*BesselJ(1,.1454...
udv := 1/Pi*((.1410265832345873e-15*BesselJ(1,.1454...
udv := 1/Pi*((.1410265832345873e-15*BesselJ(1,.1454...
udv := 1/Pi*((.1410265832345873e-15*BesselJ(1,.1454...
udv := 1/Pi*((.1410265832345873e-15*BesselJ(1,.1454...
udv := 1/Pi*((.1410265832345873e-15*BesselJ(1,.1454...
udv := 1/Pi*((.1410265832345873e-15*BesselJ(1,.1454...
udv := 1/Pi*((.1410265832345873e-15*BesselJ(1,.1454...
udv := 1/Pi*((.1410265832345873e-15*BesselJ(1,.1454...
udv := 1/Pi*((.1410265832345873e-15*BesselJ(1,.1454...
udv := 1/Pi*((.1410265832345873e-15*BesselJ(1,.1454...
udv := 1/Pi*((.1410265832345873e-15*BesselJ(1,.1454...

> plot(evalf(udv),UT=0..10); # Visibility values

[Maple Plot]

> I1 := brightRatio*r1*bessel1;

I1 := .9696273622190718e-8*BesselJ(1,.1939254724438...
I1 := .9696273622190718e-8*BesselJ(1,.1939254724438...

> I2 := r2*bessel2;

I2 := .7272205216643038e-8*BesselJ(1,.1454441043328...
I2 := .7272205216643038e-8*BesselJ(1,.1454441043328...

> u := cos(theta)*uCoord + sin(theta)*vCoord;

u := cos(37)*(16386537.89567023*sin(5.0360851640177...
u := cos(37)*(16386537.89567023*sin(5.0360851640177...

> nominatorPhase := I1*sin(Pi*u*sep) - I2*sin(Pi*u*sep);

nominatorPhase := .9696273622190718e-8*BesselJ(1,.1...
nominatorPhase := .9696273622190718e-8*BesselJ(1,.1...
nominatorPhase := .9696273622190718e-8*BesselJ(1,.1...
nominatorPhase := .9696273622190718e-8*BesselJ(1,.1...
nominatorPhase := .9696273622190718e-8*BesselJ(1,.1...
nominatorPhase := .9696273622190718e-8*BesselJ(1,.1...
nominatorPhase := .9696273622190718e-8*BesselJ(1,.1...
nominatorPhase := .9696273622190718e-8*BesselJ(1,.1...
nominatorPhase := .9696273622190718e-8*BesselJ(1,.1...

> denominatorPhase := I1*cos(Pi*u*sep) + I2*cos(Pi*u*sep);

denominatorPhase := .9696273622190718e-8*BesselJ(1,...
denominatorPhase := .9696273622190718e-8*BesselJ(1,...
denominatorPhase := .9696273622190718e-8*BesselJ(1,...
denominatorPhase := .9696273622190718e-8*BesselJ(1,...
denominatorPhase := .9696273622190718e-8*BesselJ(1,...
denominatorPhase := .9696273622190718e-8*BesselJ(1,...
denominatorPhase := .9696273622190718e-8*BesselJ(1,...
denominatorPhase := .9696273622190718e-8*BesselJ(1,...
denominatorPhase := .9696273622190718e-8*BesselJ(1,...

> phase := arctan(nominatorPhase/denominatorPhase);

phase := arctan((.9696273622190718e-8*BesselJ(1,.19...
phase := arctan((.9696273622190718e-8*BesselJ(1,.19...
phase := arctan((.9696273622190718e-8*BesselJ(1,.19...
phase := arctan((.9696273622190718e-8*BesselJ(1,.19...
phase := arctan((.9696273622190718e-8*BesselJ(1,.19...
phase := arctan((.9696273622190718e-8*BesselJ(1,.19...
phase := arctan((.9696273622190718e-8*BesselJ(1,.19...
phase := arctan((.9696273622190718e-8*BesselJ(1,.19...
phase := arctan((.9696273622190718e-8*BesselJ(1,.19...
phase := arctan((.9696273622190718e-8*BesselJ(1,.19...
phase := arctan((.9696273622190718e-8*BesselJ(1,.19...
phase := arctan((.9696273622190718e-8*BesselJ(1,.19...
phase := arctan((.9696273622190718e-8*BesselJ(1,.19...
phase := arctan((.9696273622190718e-8*BesselJ(1,.19...
phase := arctan((.9696273622190718e-8*BesselJ(1,.19...
phase := arctan((.9696273622190718e-8*BesselJ(1,.19...
phase := arctan((.9696273622190718e-8*BesselJ(1,.19...

> plot(evalf(phase), UT=0..10);

[Maple Plot]

>

>

>

>

>

>

>

>

>

>

>

>

>

>

>