maple_fringe_vis_example.html

vis_calc_ext.mws

VLTI Baseline Vector and OPD

(16.08.2007)

1. Initialisations

> restart:

> Digits := 16;

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

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

> 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)

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

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);

> # 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);

> 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)

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

> # Computes Sun coordinates

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

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

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

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

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

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

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

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

> Xs := cos(lambda0/deg2rad);

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

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

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 := LST-3*360.;

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

>

> LSTH := (LST)/15;

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

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

5. Sun altitude and azimuth

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

> HAS := LST - alpha0;

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

> ALTS := arcsin(sas);

> altsdeg := ALTS*deg2rad;

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

> # 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

>

> 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

> ds := sign(decd);

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

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

> HA := HA + 360.;

> HA/15; # Hour Angle in Hours

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

> # 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);

> altitude := arcsin(sa);

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

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

> # 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);

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

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

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

> azd := azimuth*deg2rad;

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

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

> plot(270-azdeg,UT=0..10); # VisCalc azimuth 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)

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

> # 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

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

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

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

>

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

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

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

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

> harad := HA/deg2rad;

> decrad := delta/deg2rad;

> decrad := -.2917536010865021;

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

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

> lat := latitude/deg2rad;

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

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

> 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)]]);

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

> u := UVW[1];

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

> v:= UVW[2];

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

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

>

> w := UVW[3];

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

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

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

9. Model Uniform Disk

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

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

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

> C(rho);

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

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

> evalf(uda);

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

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

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

> # 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)

> uCoord := u/lambda;

> vCoord := v/lambda;

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

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

>

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

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

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

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

10. Model Binary Disc

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

> uCoord := u/lambda;

> vCoord := v/lambda;

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

> brightRatio := 1.0; # Brightness ratio

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

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

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

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

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

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

> bessel1 := BesselJ(1,ud1);

> bessel2 := BesselJ(1,ud2);

> 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_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;