; VLTI_HRANGE - Calculate the hour angle range where an object is visible. ; hrange = vlti_hrange(declin, latitu, minalt, /silent) ; ; VLTI_BLPROJ - Calculate projected baselines positions or tracks for the VLTI. ; bltrck = vlti_blproj(declin, basein, hangin, hrange=hrange) ; ; VLTI_UVTRCK - Calculate a set of UV tracks for a certain declination. ; uvtrck = vlti_uvtrck(declin, hrs, /print) ; ; VLTI_UVCHCK - Check the UV coordinates saved in the fits header. ; result = vlti_uvchck(filnam, /silent, bltrak=bltrak) ; ; VLTI_PLTUVP - Plot a UV plane and UV tracks for a certain declination. ; vlti_pltuvp, pltset, -45, titles='UV Plane', psfile='test' ; vlti_pltuvp, pltset, /close ;+ ; NAME: ; VLTI_HRANGE ; ; PURPOSE: ; Calculate the hour angle range where an object is visible. ; ; CATEGORY: ; VLTI utilities. ; ; CALLING SEQUENCE: ; hrange = vlti_hrange(declin, latitu [, minalt] [, /silent]) ; ; REQUIRED INPUTS: ; declin declination of the target ; latitu latitude of the observatory ; ; OPTIONAL INPUTS: ; minalt minimum altitude the object must have ; ; KEYWORDS: ; silent do not print out the result ; ; OUTPUTS: ; the hour angle range where the object is higher than MINALT ; ; DESCRIPTION AND EXAMPLE: ; This programme calculates the hour angle range where the object with ; the declination DECLIN is higher than MINALT for an observatory at ; latitude LATITU. For NGC1068 at the VLT this is done by: ; hrange = vlti_hrange(-0.01344, -24.62795D, 21.0) ; ; CALLED BY: ; vlti_blproj ; vlti_pltuvt ; ; CALLING: ; None. ; ; MODIFICATION HISTORY: ; 2009-07-06 Written by Konrad R. Tristram ; FUNCTION VLTI_HRANGE, DECLIN, LATITU, MINALT, SILENT=SILENT ; SET THE MINIMUM ALTITUDE THE OBJECT MUST HAVE ;------------------------------------------------------------------------------- if n_elements(minalt) le 0 then minalt = 21.0 msgstr = 'Object never rises' + string(minalt, format='(F7.2)') + $ ' degrees over the horizon.' ; CALCULATE THE ALTITUDES FOR HOUR ANGLES BETWEEN -12 AND 0 ;------------------------------------------------------------------------------- hadec2altaz, (findgen(120)/10.-12.)/24.*360., latitu, declin, altitu ; CHECK FOR THREE DIFFERENT CASES ;------------------------------------------------------------------------------- case 1 of ; OBJECT ALWAYS OBSERVABLE: CIRCUMPOLAR OBJECT ;----------------------------------------------------------------------- min(altitu) ge minalt: hrange = [-12., +12.] ; OBJECT NEVER OBSERVABLE (NEVER RISES HIGH ENOUGH) ;----------------------------------------------------------------------- max(altitu) lt minalt: begin message, msgstr, /conti return, -1 end ; OTHER CASES: OBJECT PARTIALLY OBSERVABLE ;----------------------------------------------------------------------- else: hrange = [1, -1] # interpol(findgen(120)/10.-12., $ altitu, minalt, /spline) endcase ; PRINT OUT THE RESULT ;------------------------------------------------------------------------------- if not keyword_set(silent) then begin print, hrange, format='(" The object is observable for ' + $ 'HRANGE = [", F7.3, ",", F7.3, "].")' endif ; RETURN THE HOUR ANGLE RANGE ;------------------------------------------------------------------------------- return, hrange END ;+ ; NAME: ; VLTI_BLPROJ ; ; PURPOSE: ; Calculate projected baselines positions or tracks for the VLTI. ; ; CATEGORY: ; VLTI utilities. ; ; CALLING SEQUENCE: ; bltrck = vlti_blproj(declin, basein [, hangin] [, hrange=hrange] ; [, hdelta=hdelta] [, numele=numele]) ; ; REQUIRED INPUTS: ; declin declination of the target ; basein baseline configurations (as string array or indices) ; ; OPTIONAL INPUTS: ; hangin hour angles for which the calculation shall be done ; ; KEYWORDS: ; hrange hours before and after meridian where track begins / ends ; hdelta delta for generating uniform steps in hour angle ; numele number of elements between start and end of the UV track ; ; OUTPUTS: ; bltrck structure containing the UV tracks ; ; DESCRIPTION AND EXAMPLE: ; This programme calculates the projected baseline coordinates for an ; object at the declination DECLIN for certain baselines of the VLTI. ; The baselines can be either specified by a 2 x n integer array (every ; pair of values specifies the indices of the two stations and the total ; number of baselines is n) or by an array of strings (every element ; corresponds to a certain baseline configuration such as 'U1-U3' or ; 'E0-G0'). The baseline coordinates are calculated for the hour angles ; specified in HANGIN. If HANGIN is not provided, evenly spaced hour ; angles in the range of HRANGE are calculated either with a certain ; delta in hour angle (HDELTA) or by deviding the range into NUMELE ; steps. If HRANGE is also not provided it defaults to the time where ; the object is 21 degrees above the horizon. For example, the programme ; can be called by: ; bltrck = vlti_blproj(-45.0, ['U1-U3', 'E0-G0']) ; ; CALLED BY: ; vlti_uvtrck ; vlti_pltuvt ; ; CALLING: ; vlti_hrange ; ; MODIFICATION HISTORY: ; 2009-07-06 Written by Konrad R. Tristram ; 2011-10-21 K. Tristram: Remove fudging with PA if they decrease. ; 2012-11-17 K. Tristram: Fix calculation of PA, don't shift PAs any more. ; FUNCTION VLTI_BLPROJ, DECLIN, BASEIN, HANGIN, HRANGE=HRANGE, HDELTA=HDELTA, $ NUMELE=NUMELE ; SET A DEFAULT VALUE FOR THE NUMBER OF ELEMENTS AND MINALT ;------------------------------------------------------------------------------- if n_elements(numele) eq 0 then numele = 512 else numele = numele[0] > 2 minalt = 21. ; DEFINE THE STATION COORDINATES AND NAMES ;------------------------------------------------------------------------------- statio = [[ -16.000 , -16.000 , 0.], $ ; U1 [ 24.000 , 24.000 , 0.], $ ; U2 [ 64.0013, 47.9725, 0.], $ ; U3 [ 112.000 , 8.000 , 0.], $ ; U4 [ -32.001 , -48.013 , 0.], $ ; A0 [ -32.001 , -64.021 , 0.], $ ; A1 [ -23.991 , -48.019 , 0.], $ ; B0 [ -23.991 , -64.011 , 0.], $ ; B1 [ -23.991 , -72.011 , 0.], $ ; B2 [ -23.991 , -80.029 , 0.], $ ; B3 [ -23.991 , -88.013 , 0.], $ ; B4 [ -23.991 , -96.012 , 0.], $ ; B5 [ -16.002 , -48.013 , 0.], $ ; C0 [ -16.002 , -64.011 , 0.], $ ; C1 [ -16.002 , -72.019 , 0.], $ ; C2 [ -16.002 , -80.010 , 0.], $ ; C3 [ 0.010 , -48.012 , 0.], $ ; D0 [ 0.010 , -80.015 , 0.], $ ; D1 [ 0.010 , -96.012 , 0.], $ ; D2 [ 16.011 , -48.016 , 0.], $ ; E0 [ 32.017 , -48.017 , 0.], $ ; G0 [ 32.020 ,-112.010 , 0.], $ ; G1 [ 31.995 , -24.003 , 0.], $ ; G2 [ 64.015 , -48.007 , 0.], $ ; H0 [ 72.001 , -87.997 , 0.], $ ; I1 [ 88.016 , -71.992 , 0.], $ ; J1 [ 88.016 , -96.005 , 0.], $ ; J2 [ 88.016 , 7.996 , 0.], $ ; J3 [ 88.016 , 23.993 , 0.], $ ; J4 [ 88.016 , 47.987 , 0.], $ ; J5 [ 88.016 , 71.990 , 0.], $ ; J6 [ 96.002 , -48.006 , 0.], $ ; K0 [ 104.021 , -47.998 , 0.], $ ; L0 [ 112.013 , -48.000 , 0.]] ; M0 stanam = ['U1', 'U2', 'U3', 'U4', 'A0', 'A1', 'B0', 'B1', 'B2', 'B3', 'B4', $ 'B5', 'C0', 'C1', 'C2', 'C3', 'D0', 'D1', 'D2', 'E0', 'G0', 'G1', $ 'G2', 'H0', 'I1', 'J1', 'J2', 'J3', 'J4', 'J5', 'J6', 'K0', 'L0'] ; SET THE LATITUDE OF THE OBSERVATORY ;------------------------------------------------------------------------------- latitu = ten(-24., 37., 40.61388) ; value used by vs_baseproj latitu = -0.429838786D / !dpi *180. ; value used by EWS latitu = -24.62743941D ; value used by ESO (FITS headers) latitu = -24.62794830D ; value used by ESO (web site) colati = (90. + latitu) * !dpi / 180. ; latitude from the South Pole ; IF NO HOUR ANGLES ARE PROVIDED GENERATE THEM HERE ;------------------------------------------------------------------------------- if n_elements(hangin) le 0 then begin ; IF HRANGE NOT GIVEN, CALCULATE HOUR ANGLE WHERE OBJECT IS OBSERVABLE ;----------------------------------------------------------------------- if n_elements(hrange) ne 2 then begin hrange = vlti_hrange(declin, latitu, minalt) if hrange[0] eq -1 then return, -1 endif else hrange = float(hrange) ; IF HDELTA GIVEN, CALCULATE HOUR ANGLES SPACED BY HDELTA ;----------------------------------------------------------------------- if n_elements(hdelta) ge 1 then begin tmpvar = findgen(fix(hrange[1]/hdelta[0]) + 1) * hdelta[0] if n_elements(tmpvar) eq 1 then hangle = tmpvar else $ hangle = [-reverse(tmpvar), tmpvar[1:*]] ; OTHERWISE SPACE HOUR ANGLES UNIFORMLY ;----------------------------------------------------------------------- endif else hangle = hrange[0] + (hrange[1]-hrange[0]) * $ (findgen(numele)/(numele-1)) endif else hangle = hangin ; DETERMINE HOW MANY HOUR ANGLES THERE ARE FINALLY ;------------------------------------------------------------------------------- numele = n_elements(hangle) ; IF BASEIN IS A TWO DIMENSIONAL ARRAY OF NUMBERS THEN USE IT AS STATION INDICES ;------------------------------------------------------------------------------- if total(size(basein, /type) eq [1,2,3,4,5]) and $ ((size(basein))[0] eq 2) then begin baseli = fix(basein) ; IF BASEIN IS A STRING ARRAY THEN INTERPRET IT AS BASELINE NAMES ;------------------------------------------------------------------------------- endif else if size(basein, /type) eq 7 then begin baseli = [0] ; LOOP OVER ALL ELEMENTS IN THE STRING ARRAY ;----------------------------------------------------------------------- for i=0,n_elements(basein)-1 do begin ; SPLIT THE STRING INTO TWO PARTS ;--------------------------------------------------------------- tmpvar = strsplit(basein[i], '-', /extract) ; MATCH THE POSITION STRINGS TO THE STATION NAMES ;--------------------------------------------------------------- tmpvar = [where(tmpvar[0] eq stanam), $ where(tmpvar[1] eq stanam)] ; COPY STATION INDICES ONLY IF A VALID BASELINE WAS FOUND ;--------------------------------------------------------------- if (tmpvar[0] ge 0) and (tmpvar[1] ge 0) then begin baseli = [baseli, tmpvar] endif else message, 'Skipping invalid baseline string.', /contin endfor ; EXIT IF NO VALID BASELINE NAMES WHERE FOUND ;----------------------------------------------------------------------- if n_elements(baseli) le 1 then begin message, 'No valid baseline station array provided.', /continue return, -1 endif else baseli = reform(baseli[1:*], 2, (n_elements(baseli)-1)/2) ; OTHERWISE DON'T KNOW HOW TO INTERPRET BASEIN, THEREFORE ALSO EXIT ;------------------------------------------------------------------------------- endif else begin message, 'No valid baseline station array provided.', /continue return, -1 endelse ; CREATE THE STRUCT HOLDING THE RESULT ;------------------------------------------------------------------------------- numbli = n_elements(baseli[0,*]) bltrck = {baseli:'', hangle:hangle, $ ucoord:fltarr(numele), vcoord:fltarr(numele), $ pbline:fltarr(numele), pangle:fltarr(numele)} bltrck = replicate(bltrck, numbli) ; ROTATE THE VLTI COORDINATES INTO GEOGRAPHICAL COORDINATES ;------------------------------------------------------------------------------- vltrot = 18.984 * !dpi / 180. vltrot = [[cos(vltrot), -sin(vltrot)], [sin(vltrot), cos(vltrot)]] for i=0,n_elements(statio[0,*])-1 do statio[0:1, i] = vltrot ## statio[0:1, i] ; CALCULATE THE BASELINE COORDINATES ;------------------------------------------------------------------------------- blcoor = fltarr(3, numbli) blcoor[0, *] = statio[0, baseli[1, *]] - statio[0, baseli[0, *]] ; east blcoor[1, *] = statio[1, baseli[1, *]] - statio[1, baseli[0, *]] ; north blcoor[2, *] = statio[2, baseli[1, *]] - statio[2, baseli[0, *]] ; elevation ; CHANGE HOUR ANGLE INTO RADIAN INSTEAD OF HOURS ;------------------------------------------------------------------------------- hangle = 15. * !dpi / 180. * hangle ; LOOP OVER THE BASELINES ;------------------------------------------------------------------------------- for i=0, numbli-1 do begin ; CALCULATE PROJECTIONS NORTH, WEST AND IN ALTITUDE ;----------------------------------------------------------------------- projen = blcoor[1, i]*cos(colati) + blcoor[2, i]*sin(colati) ; north projew = -blcoor[0, i] ; west projea = -blcoor[1, i]*sin(colati) + blcoor[2, i]*cos(colati) ; altitude ; ADD HOUR ANGLE DEPENDENCE AND COPY TO THE RESULT STRUCT ;----------------------------------------------------------------------- bltrck[i].ucoord = projen*sin(hangle) - projew*cos(hangle) bltrck[i].vcoord = -projen*sin(declin*!dpi/180.)*cos(hangle) - $ projew*sin(declin*!dpi/180.)*sin(hangle) - $ projea*cos(declin*!dpi/180.) ; CALCULATE THE BASELINE LENGTH AND POSITION ANGLE ;----------------------------------------------------------------------- bltrck[i].pbline = sqrt(bltrck[i].ucoord^2 + bltrck[i].vcoord^2) bltrck[i].pangle = atan(bltrck[i].ucoord, bltrck[i].vcoord) / !dpi*180. ; MOVE THE POSITION ANGLES TO THE FIRST TWO QUADRANTS ;----------------------------------------------------------------------- ; This should not be done if we also consider the phases, as they ; will switch sign if we move the position angles or interchange the ; telescopes!!! ; sel_pa = where(bltrck[i].pangle lt -60., count) ; if (count gt 0) then begin ; bltrck[i].pangle[sel_pa] = bltrck[i].pangle[sel_pa] + 180. ; endif ; COPY THE BASELINE NAME TO THE RESULT STRUCT ;----------------------------------------------------------------------- bltrck[i].baseli = stanam[baseli[0,i]]+'-'+stanam[baseli[1,i]] endfor ; RETURN THE RESULT ;------------------------------------------------------------------------------- return, bltrck END ;+ ; NAME: ; VLTI_UVTRCK ; ; PURPOSE: ; Calculate a set of UV tracks of the VLTI for a certain declination. ; ; CATEGORY: ; VLTI utilities. ; ; CALLING SEQUENCE: ; uvtrck = vlti_uvtrck(declin, hrs [, /print]) ; ; REQUIRED INPUTS: ; dec declination of the target ; hrs hours before and after meridian where track begins / ends ; ; OPTIONAL INPUTS: ; none ; ; KEYWORDS: ; print print out the uv tracks into an ASCII file ; ; OUTPUTS: ; uvtrck structure containing the uv tracks ; ; DESCRIPTION AND EXAMPLE: ; This programme calculates a set of UV tracks of the VLTI for a ; certain declination. The tracks are calculated symmetrically to the ; culmination of the object with the number of hours on either side ; specified in HRS. The programme is invoked by e.g.: ; uvtrck = vlti_uvtrck( -65.33, 6.0) ; ; CALLED BY: ; none ; ; CALLING: ; vlti_blproj ; (vs_baseproj) ; ; MODIFICATION HISTORY: ; 2007-04-27 Written by Konrad R. Tristram ; 2009-07-06 K. Tristram: Renamed and rewritten using VLTI_BLPROJ. ; 2012-11-17 K. Tristram: Fix calculation of PA, don't shift PAs any more. ; FUNCTION VLTI_UVTRCK, DECLIN, HRS, PRINT=PRINT ; SET NUMBER OF ELEMENTS IN THE TRACK ;------------------------------------------------------------------------------- numelm = 512 ; INITIALISE THE BASELINES TO BE CALCULATED ;------------------------------------------------------------------------------- blname = ['U1-U2', 'U2-U3', 'U3-U4', 'U1-U3', 'U2-U4', 'U1-U4', 'E0-G0', $ 'G0-H0', 'A0-B1', 'A0-B2', 'A0-C1', 'A1-B0', 'A1-B2', 'A1-C0', $ 'A1-C2', 'B0-D0', 'B2-C0', 'B2-C1'] baseli = [ [0,1], [1,2], [2,3], [0,2], [1,3], [0,3], [19,20], $ [20,23], [ 4, 7], [ 4, 8], [ 4,13], [ 5, 6], [ 5, 8], [ 5,12], $ [ 5,14], [ 6,16], [ 8,12], [ 8,13]] ; CALCULATE THE BASELINE TRACKS ;------------------------------------------------------------------------------- uvtrck = vlti_blproj(declin, blname, hrange=[-1.0, +1.0]*hrs[0]) ; ALTERNATIVELY USE VS_BASEPROJ TO CALCULATE THE UV TRACKS (OLD METHOD) ;------------------------------------------------------------------------------- if 0 then begin ; COMPILE THE FUNCTION VS_BASEPROJ ;----------------------------------------------------------------------- forward_function vs_baseproj cd, '~/progs/VS/', curren=curren resolve_routine, "vs_baseproj", /compile_full_file, /is_function cd, curren ; RUN THE FUNCTION VS_BASEPROJ ;----------------------------------------------------------------------- hangle = findgen(numelm)/numelm*2.*hrs - hrs bltrak = vs_baseproj(declin, -1.0 * hrs, +1.0 * hrs, baseli, $ n_elements(blname), 4.84814e-6, numelm) ; CREATE THE STRUCT THAT WILL HOLD THE RESULT ;----------------------------------------------------------------------- nularr = make_array(numelm, /float, value=0.0) hangle = findgen(numelm)/numelm*2.*hrs - hrs uvtrck = {uvtrck, baseli:'', hangle:hangle, ucoord:nularr, $ vcoord:nularr, pbline:nularr, pangle:nularr} uvtrck = replicate(uvtrck, n_elements(blname)) ; COPY THE DATA INTO THE RESULT ;----------------------------------------------------------------------- uvtrck.ucoord = transpose(-bltrak[*,*,0]) uvtrck.vcoord = transpose(+bltrak[*,*,1]) uvtrck.pbline = sqrt(uvtrck.ucoord^2 + uvtrck.vcoord^2) uvtrck.pangle = atan(uvtrck.ucoord, uvtrck.vcoord)/!dpi*180. endif ; PRINT OUT THE RESULT TO A FILE IF KEYWORD IS SET ;------------------------------------------------------------------------------- if keyword_set(print) then begin for i=0,n_elements(blname)-1 do begin openw, 1, blname[i] + '.dat' for j=0,numelm-1 do begin printf, 1, uvtrck[i].ucoord[j], uvtrck[i].vcoord[j] endfor free_lun, 1 endfor endif ; RETURN THE RESULT ;------------------------------------------------------------------------------- return, uvtrck END ;+ ; NAME: ; VLTI_UVCHCK ; ; PURPOSE: ; Check the UV coordinates saved in the fits header. ; ; CATEGORY: ; VLTI utilities. ; ; CALLING SEQUENCE: ; result = vlti_uvchck(filnam [, /silent] [, bltrak=bltrak]) ; ; REQUIRED INPUTS: ; filnam name of a file for which to check the UV coordinates ; ; OPTIONAL INPUTS: ; none ; ; KEYWORDS: ; silent do not print out the result ; bltrak structure containing the recalculated baseline data ; ; OUTPUTS: ; result 1 if baseline data is within errors, 0 otherwise ; ; DESCRIPTION AND EXAMPLE: ; This function checks the UV coordinates stored in the FITS header of ; an interferometry file by recalculating the projected baseline length ; and the position angle. Optionally the correct baseline data for the ; start and the end of the data recording (fringe track) can be obtained ; through the keyword BLTRAK. The function currently only works for MIDI ; FITS files. The programme is used by e.g.: ; result = vlti_uvchck('MIDI.2009-01-01T12:34:56.000.fits') ; ; CALLED BY: ; midi_getdat ; ; CALLING: ; vlti_blproj ; ; MODIFICATION HISTORY: ; 2007-04-27 Written by Konrad R. Tristram ; 2009-07-06 K. Tristram: Rewritten completely using VLTI_BLPROJ. ; 2011-10-21 K. Tristram: Fixed calculation of PBLINE for EWS. ; 2012-11-17 K. Tristram: Fix calculation of PANGLE for EWS. ; FUNCTION VLTI_UVCHCK, FILNAM, SILENT=SILENT, BLTRAK=BLTRAK ; READ IN THE FITS HEADER OF THE FILE ;------------------------------------------------------------------------------- filarr = strsplit(filnam[0], /extr) numfil = n_elements(filarr) dummyv = readfits(filarr[0], header, /silent) ; READ THE FITS HEADER AND GET HEADER VALUES NEEDED FOR THE CALCULATION ;------------------------------------------------------------------------------- objnam = oirGetKey(filarr[0], 'OBS TARG NAME') rascen = sxpar(header, 'RA') declin = sxpar(header, 'DEC') exptim = sxpar(header, 'EXPTIME')/3600. lstfil = sxpar(header, 'LST')/3600. + [0, exptim] ; DETERMINE PRECISE LST AT START AND END OF THE OBSERVATION ;------------------------------------------------------------------------------- if ~ keyword_set(quick) then begin ; GET THE GEOGRAPHIC LONGITUDE OF THE OBSERVATORY ;----------------------------------------------------------------------- geolon = oirGetKey(filarr[0], 'ISS GEOLON') ; READ THE ALL DATA ;----------------------------------------------------------------------- dummyv = oirgetdata(filnam) ; CONVERT JULIAN DATE TO CIVIL DATE AND THEN TO LST ;----------------------------------------------------------------------- daycnv, reform(dummyv[[0,n_elements(dummyv)-1]].time + 2400000.5D), $ dat_yr, dat_mn, dat_dy, dat_hr ct2lst, lstime, geolon, 0, dat_hr, dat_dy, dat_mn, dat_yr endif else lstime = lstfil ; CALCULATE THE HOUR ANGLE ;------------------------------------------------------------------------------- hangle = reform(lstime) - rascen/360.*24. ; GET THE BASELINE NAME ;------------------------------------------------------------------------------- blname = strcompress(oirGetKey(filarr[0], 'CONF STATION1'), /remove) + '-' + $ strcompress(oirGetKey(filarr[0], 'CONF STATION2'), /remove) ; CALCULATE THE UV COORDINATES FOR THE TIME OF OBSERVATIONS WITH VLTI_BLPROJ ;------------------------------------------------------------------------------- bltrak = vlti_blproj(declin, blname, hangle) our_uv = [[bltrak.ucoord], [bltrak.vcoord]] our_bl = [[bltrak.pbline], [bltrak.pangle]] ; USE EWS TO CALCULATE THE BASELINE COORDINATES ;------------------------------------------------------------------------------- ews_uv = midifitsuvw(filarr[0]) ews_bl = [sqrt(total(ews_uv[0:1]^2)), atan(ews_uv[0], ews_uv[1])/!dpi*180.] ; GET THE PROJECTED BASELINE LENGTH AND POSITION ANGLE ESO CALCULATED ;------------------------------------------------------------------------------- e_pbli = [oirGetKey(filarr[numfil-1], 'PBL12 START'), $ oirGetKey(filarr[numfil-1], 'PBL12 END')] e_pang = [oirGetKey(filarr[numfil-1], 'PBLA12 START'), $ oirGetKey(filarr[numfil-1], 'PBLA12 END')] ; REMOVE 180 DEGREE AMBIGUITY IN THE POSITION ANGLE ;------------------------------------------------------------------------------- fixidx = where((e_pang - our_bl[*,1]) gt 90.) if fixidx[0] ge 0 then begin our_bl[fixidx,1] += 180. endif fixidx = where((e_pang - our_bl[*,1]) lt -90.) if fixidx[0] ge 0 then begin our_bl[fixidx,1] -= 180. endif ; DETERMINE THE RESULT ;------------------------------------------------------------------------------- dif_pb = total(abs(our_bl[*,0]-e_pbli)) / 2. dif_pa = total(abs(our_bl[*,1]-e_pang)) / 2. result = (dif_pb le 1.0) and (dif_pa le 1.0) ; PRINT THE RESULTS ;------------------------------------------------------------------------------- if not keyword_set(silent) then begin print, 'Checking UV coordinates for ' + filarr[0] + ':' print, objnam , format='(" Object = ", A14)' print, rascen, declin, format='(" RA = ", F10.5, ", DEC = ", F10.5)' print, '====================================================' print, ' Start End' print, lstfil , format='("LST (FILE) :", 2F10.3)' print, lstime , format='("LST (FRAMES) :", 2F10.3)' print, hangle , format='("HA (FRAMES) :", 2F10.3)' print, ' print, our_uv[*,0], format='("UCOORD (CALC):", 2F10.3)' print, ews_uv[0] , format='("UCOORD (EWS) :", 1F10.3)' print, ' print, our_uv[*,1], format='("VCOORD (CALC):", 2F10.3)' print, ews_uv[1] , format='("VCOORD (EWS) :", 1F10.3)' print, ' print, e_pbli , format='("PBLINE (ESO) :", 2F10.3)' print, our_bl[*,0], format='("PBLINE (CALC):", 2F10.3)' print, ews_bl[0] , format='("PBLINE (EWS) :", 2F10.3)' print, '' print, e_pang , format='("PANGLE (ESO) :", 2F10.3)' print, our_bl[*,1], format='("PANGLE (CALC):", 2F10.3)' print, ews_bl[1] , format='("PANGLE (EWS) :", 2F10.3)' print, '' if result then print, '--> Deviations are within the errors.' $ else print, '--> Deviations significant! Header erroneous?' endif ; RETURN THE RESULT ;------------------------------------------------------------------------------- return, result END ;+ ; NAME: ; VLTI_PLTUVP ; ; PURPOSE: ; Plot a UV plane and plot UV tracks for a certain declination. ; ; CATEGORY: ; VLTI utilities. ; ; CALLING SEQUENCE: ; vlti_pltuvp, pltset [, declin] [, blname=blname] [, minalt=minalt] ; [, pblran=pblran] [, titles=titles] [, pltpos=pltpos] ; [, barpos=barpos] [, /noline] [, /noarrow] [, /nolabel] ; [, /oldwin] [, psfile=psfile] [, /close] ; ; REQUIRED INPUTS: ; pltset structure containing old and new device settings ; ; OPTIONAL INPUTS: ; declin declination of the target ; ; KEYWORDS: ; blname baseline configurations (as string array or indices) ; minalt minimum altitude the object must have ; pblran range of baseline lengths ; titles title string for the plot ; noline don't plot the cross hairs ; noarrow don't plot the arrows ; nolabel don't label the UV tracks ; oldwin reuse the old plotting window (only for screen output) ; pltsiz size of the plot (A&A single-column, 10cm, etc.) ; psfile name of the PS file for PS output ; close close the currently active plot ; ; OUTPUTS: ; pltpos position of the plot in device coordinates ; barpos position of the colour bar in device coordinates ; ; DESCRIPTION AND EXAMPLE: ; This programme sets up a plot of the UV plane. First the axes are ; plotted according to the ranges specified in PBLRAN, then a cross ; hair and arrows indicating the North and East direction a plotted ; if not suppressed by the respective keywords. If a declination is ; specified as a second parameter, UV tracks for this declination are ; overplotted in grey. First the UV tracks are calculated for the ; baselines specified in BLNAME. If BLNAME is not set then the four UT ; baselines are plotted. The UV tracks are calculated for at least an ; elevation of 21° above the horizon, if not specified differently by ; the keyword MINALT. Labels are attached at the start of the tracks ; unless suppressed by /NOLABEL. To restore the original plot settings ; and properly close the plot, the programme has to be called again ; using the keyword /CLOSE. ; Example of plotting on UV point: ; vlti_pltuvp, pltset, -45, title='UV plane' ; oplot, [50,-50], [30,-30], psym=4, color=250 ; vlti_pltuvp, pltset, /close ; ; CALLED BY: ; midi_pltuvp ; ; CALLING: ; vlti_hrange ; vlti_blproj ; ; MODIFICATION HISTORY: ; 2010-12-03 Written by Konrad R. Tristram ; 2012-03-15 K. Tristram: Suppress labels if outside the plot area. ; 2012-03-29 K. Tristram: Added keyword PLTSIZ to set size of the plot. ; 2012-11-01 K. Tristram: Adjust PS font size for !p.font=0. ; PRO VLTI_PLTUVP, PLTSET, DECLIN, BLNAME=BLNAME, MINALT=MINALT, PBLRAN=PBLRAN, $ TITLES=TITLES, PLTPOS=PLTPOS, BARPOS=BARPOS, NOLINE=NOLINE, $ NOARROW=NOARROW, NOLABEL=NOLABEL, OLDWIN=OLDWIN, $ PLTSIZ=PLTSIZ, PSFILE=PSFILE, CLOSE=CLOSE ; SET DEFAULT VALUES FOR THE BASELINE RANGE ;------------------------------------------------------------------------------- if n_elements(pblran) ne 2 then pltran=[-100.,+100.] else pltran=float(pblran) ; SET DEFAULT VALUES FOR THE OBSERVATORY LATITUDE AND BASELINES ;------------------------------------------------------------------------------- if n_elements(latitu) le 0 then latitu = -24.62794830D ; ESO's VLTI coordinates if n_elements(blname) le 0 then blname = ['U1-U2', 'U2-U3', 'U3-U4', 'U1-U3'] ; CHECK IF PLTSET WILL BE RETURNED; OTHERWISE AUTOMATICALLY CLOSE THE PLOT ;------------------------------------------------------------------------------- if ~ arg_present(pltset) then autoclose = 1B else autoclose = 0B ; CHECK IF PROGRAMME IS CALLED FOR CLOSING A PLOT ;------------------------------------------------------------------------------- if keyword_set(close) then begin ; CLOSE THE PLOTTING DEVICE AND RESET THE PLOT PARAMETERS ;----------------------------------------------------------------------- if ((!d.name eq 'PS') or (!d.name eq 'METAFILE')) then device, /close set_plot, pltset.d !p = pltset.p !x = pltset.x !y = pltset.y tmpvar = temporary(pltset) ; OTHERWISE OPEN A NEW PLOT ;------------------------------------------------------------------------------- endif else begin ; IF PLTSET IS ALREADY CORRECTLY DEFINED THEN CLOSE THE PLOT STILL OPEN ;----------------------------------------------------------------------- if size(pltset, /type) eq 8 then begin if array_equal(tag_names(pltset), ['D', 'P']) then begin message, 'A plot is still seems to be open! ' + $ 'Closing it first.', /informa, /continue vlti_pltuvp, pltset, /close endif endif ; SAVE THE OLD PLOT SETTINGS IN THE PLTSET STRUCTURE ;----------------------------------------------------------------------- pltset = {d:!d.name, p:!p, x:!x, y:!y} ; SET DEVICE PROPERTIES AND PLOT PROPERTIES FOR PS OUTPUT ;----------------------------------------------------------------------- if keyword_set(psfile) then begin set_plot, 'PS' if arg_present(barpos) then resize=1000 else resize=0 if n_elements (pltsiz) lt 1 then pltsiz=0 case pltsiz of ; DEVICE PROPERTIES FOR AN A&A 12CM WIDE PLOT ;--------------------------------------------------------------- 0: begin device, file=psfile+'.ps', xoffset= 4.5, yoffset= 9.0, $ xsize=12., ysize=12.-0.001*resize, /portrait, $ encapsulated=encaps, color=1 pltpos = [1400, 900, 1400+10200-resize,900+10200-resize] end ; DEVICE PROPERTIES FOR A 10CM WIDE PLOT ;--------------------------------------------------------------- 1: begin device, file=psfile+'.ps', xoffset= 5.5, yoffset=10.0, $ xsize=10., ysize=10.-0.001*resize, /portrait, $ encapsulated=encaps, color=1 pltpos = [1400, 900, 1400+8200-resize,900+8200-resize] end ; DEVICE PROPERTIES FOR AN A&A 8CM WIDE PLOT ;--------------------------------------------------------------- else: begin device, file=psfile+'.ps', xoffset= 6.5, yoffset=10.4, $ xsize=8.8, ysize=8.8-0.001*resize, /portrait, $ encapsulated=encaps, color=1 pltpos = [1350, 900, 1350+7100-resize, 900+7100-resize] end endcase barpos = pltpos[[2,1,2,3]] + [ 0, 0, 300, 0] if !p.font eq 0 then !p.charsize=0.65 else !p.charsize=0.75 ; SET DEVICE PROPERTIES AND PLOT PROPERTIES FOR METAFILE OUTPUT ;----------------------------------------------------------------------- endif else if keyword_set(metafile) then begin set_plot, 'METAFILE' device, file=metafile+'.emf', set_font='Arial*18', $ xsize=36.0, ysize=36.0 !p.font=0 !p.charsize=0.80 offset= 255 factor = -1 pltpos = 38*[ 3.5, 2.0, 35.5, 34.0] ; OTHERWISE SET DEVICE PROPERTIES AND PLOT PROPERTIES FOR SCREEN OUTPUT ;----------------------------------------------------------------------- endif else begin if !version.os_family eq 'Windows' then set_plot, 'WIN' $ else set_plot, 'X' if keyword_set(oldwin) then newwin=0 else newwin=1 if arg_present(barpos) then resize=50 else resize=0 window, xsize=800.0, ysize=800.0-resize, free=newwin device, decompose=0 !p.charsize=1.25 pltpos = [ 75, 50, 75+700-resize, 50+700-resize] barpos = pltpos[[2,1,2,3]] + [ 0, 0, 20, 0] endelse ; CREATE THE PLOT ;----------------------------------------------------------------------- loadct, 39, /silent plot, -pltran, pltran, xrange=-pltran, xstyle=1, ystyle=1, $ xtitle='U [m]', ytitle='V [m]', title=titles, position=pltpos, $ /device, /nodata ; PLOT HELPER LINES IF NOT SUPPRESSED ;----------------------------------------------------------------------- if ~ keyword_set(noline) then oplot, [-1D9,+1D9], [ 0.0, 0.0], linest=2 if ~ keyword_set(noline) then oplot, [ 0.0, 0.0], [-1D9,+1D9], linest=2 ; PLOT ARROWS IF NOT SUPPRESSED ;----------------------------------------------------------------------- if ~ keyword_set(noarrow) then begin tmpvar = pltpos[2] - pltpos[0] tmpvar = [pltpos[0] + 0.15*tmpvar, $ pltpos[1] + 0.85*tmpvar, 0.07*tmpvar] arrow, tmpvar[0], tmpvar[1], tmpvar[0], tmpvar[1]+tmpvar[2] arrow, tmpvar[0], tmpvar[1], tmpvar[0]-tmpvar[2], tmpvar[1] xyouts, tmpvar[0], tmpvar[1]+1.15*tmpvar[2], 'N', align=0.5, $ charsize=1.25*!p.charsize, /device xyouts, tmpvar[0]-1.15*tmpvar[2], tmpvar[1]-0.08*tmpvar[2], $ 'E', align=1.0, charsize=1.25*!p.charsize, /device endif ; CALCULATE THE UV TRACKS ;----------------------------------------------------------------------- if n_elements(declin) gt 0 then begin hrange = vlti_hrange(declin, latitu, minalt, /silent) endif else hrange = -1 if (hrange[0] ne -1) then begin $ bltrck = vlti_blproj(declin, blname, hrange=hrange, numele=256) &$ endif ; PLOT THE UVTRACKS ;----------------------------------------------------------------------- loadct, 0, /silent colour = 210 for i=0,n_elements(bltrck)-1 do begin oplot, +bltrck[i].ucoord, +bltrck[i].vcoord, color=colour oplot, -bltrck[i].ucoord, -bltrck[i].vcoord, color=colour endfor ; PLOT THE NAMES OF THE TRACKS ;----------------------------------------------------------------------- if ~ keyword_set(nolabel) then for i=0,n_elements(bltrck)-1 do begin pangle = atan((bltrck[i].ucoord[0]-bltrck[i].ucoord[1]), $ (bltrck[i].vcoord[0]-bltrck[i].vcoord[1])) if pangle lt 0 then pangle += !pi if keyword_set(psfile) then chrsiz = 0.50*!p.charsize $ else chrsiz = 0.75*!p.charsize if (abs(bltrck[i].ucoord[0]) gt pltran[0]) and $ (abs(bltrck[i].ucoord[0]) lt pltran[1]) and $ (abs(bltrck[i].vcoord[0]) gt pltran[0]) and $ (abs(bltrck[i].vcoord[0]) lt pltran[1]) then $ xyouts, bltrck[i].ucoord[0], bltrck[i].vcoord[0], $ bltrck[i].baseli, color=colour, align=0.5, $ charsize=0.75*!p.charsize, $ orient=pangle/!pi*180.-90. endfor loadct, 39, /silent ; AUTOMATICALLY CLOSE THE PLOT IF NO PLTSET WAS GIVEN ;----------------------------------------------------------------------- if autoclose then vlti_pltuvp, pltset, /close endelse END