%%%%%%%%%%%%%%%5 Function plotOvroSurf(d,channel) takes the data stored in the structure d, coordinate transform from antenna az,el to apparent az,el to apparent right ascension and declination Then grids the data in 0.2 degree bins and mesh plots the data with color scale returns the grid data xi,yi,zi that can be used with the matlab mesh function. channel is a number in range 1:6, which selects the column of data to map CJC-9 April 2010 %%%%%%%%%%%%%%%%%%%%% ogk edit 1: added channel option %%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
0001 function [xi,yi,zi] = plotOvroSurf(d,channel) 0002 %%%%%%%%%%%%%%%%5 0003 %Function plotOvroSurf(d,channel) 0004 %takes the data stored in the structure d, coordinate transform from 0005 %antenna az,el to apparent az,el to apparent right ascension and 0006 %declination 0007 %Then grids the data in 0.2 degree bins and mesh plots the data with color 0008 %scale 0009 %returns the grid data xi,yi,zi that can be used with the matlab mesh 0010 %function. 0011 % channel is a number in range 1:6, which selects the column of data to map 0012 %CJC-9 April 2010 0013 % 0014 %%%%%%%%%%%%%%%%%%%%%% 0015 % ogk edit 1: 0016 % added channel option 0017 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5 0018 0019 0020 long=-118.2822; 0021 lat=37.2339; 0022 %get the antenna coordinate first 0023 azprime = d.antenna0.servo.az; 0024 elprime = d.antenna0.servo.el; 0025 jd=mjd2jd(d.antenna0.receiver.utc); 0026 lst1=lst(jd,deg2rad(long),'m'); 0027 lst1=lst1*360; 0028 %use the iterative approach to get rid of the refraction and pointing model 0029 %and arrive back at apparent coordinates. 0030 [az,el,daz,del]=apparentAzEl_ovro_surf(d); 0031 [equa] = horiz_coo_ovro_surf([deg2rad(az) deg2rad(el)],jd,[deg2rad(long) deg2rad(lat)],'e'); 0032 equa=rad2deg(equa); 0033 %define variables 0034 0035 0036 0037 %save the data to a matrix to remove the nan's 0038 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0039 % ogk edit 1: 0040 % added channel option here 0041 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0042 dmap2=[equa d.antenna0.receiver.data(:,channel)]; 0043 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0044 0045 %remove the nans 0046 dmap2=cleanRaDec_ovro_surf(dmap2); 0047 0048 0049 [rai,deci]=meshgrid(min(dmap2(:,1)):0.2:max(dmap2(:,1)),min(dmap2(:,2)):0.2:max(dmap2(:,2))); 0050 %deci=56:0.3:59.5; 0051 [xi,yi,zi] = griddata(dmap2(:,1),dmap2(:,2),-dmap2(:,3),rai,deci,'nearest'); 0052 mesh(xi,yi,zi); 0053 title('Grid Data of Scan'); 0054 xlabel('Right Ascension [deg]'); 0055 ylabel('Declination [deg]'); 0056 colorbar 0057 axis square 0058 0059 0060 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5 0061 function OutCoo=horiz_coo_ovro_surf(InCoo,JD,TopoPos,Direction); 0062 %-------------------------------------------------------------------- 0063 % horiz_coo function Horizontal coordinates conversion 0064 % Converting from equatorial coordinates 0065 % to Horizontal coordinates and visa 0066 % versa. 0067 % input : - Two columns matrix of coordinates, 0068 % (Long & Lat) | (Az & Alt) in radians 0069 % - vector of JDs + UT fraction of day, 0070 % if scalar value is given, then it duplicated 0071 % for all coordinates. 0072 % - Geodetic Coordinates, east long & north lat in radians 0073 % if scalar value is given then it is duplicate for 0074 % all coordinates. 0075 % - Direction, 0076 % 'h' - from equatorial to horizontal (default). 0077 % 'e' - from horizontal to equatorial. 0078 % output : - two column matrix of output coordinates. 0079 % Tested : matlab 5.3 0080 % By : Eran O. Ofek August 1999 0081 % URL : http://wise-obs.tau.ac.il/~eran/matlab.html 0082 %-------------------------------------------------------------------- 0083 if (nargin==3), 0084 Direction = 'h'; 0085 elseif (nargin==4), 0086 % no default 0087 else 0088 error('Illigal number of input arguments'); 0089 end 0090 0091 N = length(InCoo(:,1)); 0092 if (length(JD)==1), 0093 JD = JD.*ones(N,1); 0094 end 0095 0096 if (length(TopoPos(:,1))==1), 0097 TopoPos = ones(N,1)*TopoPos; 0098 end 0099 0100 0101 % Don't convert Geodetic latitude to Geocentric. 0102 %GeodLat = TopoPos(:,2); 0103 %GeocLatTemp = geod2geoc([TopoPos(:,1),GeodLat,zeros(N,1)],'WGS84'); 0104 %GeocLat = GeocLatTemp(:,2); 0105 %TopoPos(:,2) = GeocLat; 0106 0107 0108 % calculating Local Mean Sidereal Time 0109 LST=lst(JD,TopoPos(:,1),'m'); 0110 0111 0112 if (Direction=='h'), 0113 % convert equatorial to horizontal 0114 0115 % calculate the Hour Angle 0116 HA = LST.*2.*pi - InCoo(:,1); 0117 Dec = InCoo(:,2); 0118 Lat = TopoPos(:,2); 0119 0120 SinAlt = sin(Dec).*sin(Lat) + cos(Dec).*cos(HA).*cos(Lat); 0121 CosAlt = sqrt(1-SinAlt.*SinAlt); 0122 0123 SinAz = (-cos(Dec).*sin(HA))./CosAlt; 0124 CosAz = (sin(Dec).*cos(Lat) - cos(Dec).*cos(HA).*sin(Lat))./CosAlt; 0125 0126 Az = atan2(SinAz, CosAz); 0127 Alt = asin(SinAlt); 0128 0129 I = find(Az<0); 0130 Az(I) = 2.*pi+Az(I); 0131 0132 OutCoo = [Az, Alt]; 0133 elseif (Direction=='e'), 0134 Az = InCoo(:,1); 0135 Alt = InCoo(:,2); 0136 Lat = TopoPos(:,2); 0137 0138 SinDec = sin(Alt).*sin(Lat) + cos(Alt).*cos(Az).*cos(Lat); 0139 CosDec = sqrt(1 - SinDec.*SinDec); 0140 0141 SinHA = (-cos(Alt).*sin(Az))./CosDec; 0142 CosHA = (sin(Alt).*cos(Lat) - cos(Alt).*cos(Az).*sin(Lat))./CosDec; 0143 HA = atan2(SinHA, CosHA); 0144 0145 RA = LST.*2.*pi - HA; 0146 Dec = asin(SinDec); 0147 0148 % converting to range [0,1) 0149 RA = 2.*pi.*(RA./(2.*pi) - floor(RA./(2.*pi))); 0150 0151 OutCoo = [RA, Dec]; 0152 else 0153 error('Illigal type of conversion'); 0154 end 0155 0156 0157 0158 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0159 function [aza,ela,erraz,errel] = apparentAzEl_ovro_surf(d) 0160 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0161 %very rough function to go from antenna coordinates back to apparent az/el 0162 %in preparation to go to RA/DEC-CJC 0163 %d is the C-BASS data structure 0164 %aza is the apparent azimuth 0165 %ela is the apparent elevation 0166 %erraz is the error between the antenna coordinate recorder in the data and 0167 %what we get when we use the apparent position to start with for refraction 0168 %and pointing calcualtions 0169 %similarly for erralt 0170 %CJC 0171 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5 0172 0173 0174 0175 %get the antenna coordinates that we start with 0176 azprime=d.antenna0.servo.az; 0177 elprime = d.antenna0.servo.el; 0178 0179 %lets just see what sort of correction we get at this antenna angle (AZ')- 0180 %if we assume the difference between AZ (apparent) and AZ'(antenna) is 0181 %small then we can use this number as a start approximation- here we have 3 0182 %coordinates Az (apparent AZ) Az' (antenna Az) AZ'' (meaningless azimuth 0183 %position calculated by applying a pointing and refraction correction to Az') and we 0184 %make the assumptio that Az-Az' is approximately Az'-Az'' 0185 [azdoubleprime,eldoubleprime]=calcAzEl2(azprime,elprime,d); 0186 %now we calculate what sort of offset we have, daz and del-note that this 0187 %offset is between AZ' and AZ''- ie daz=AZ'-AZ''- therefore to go back to 0188 %AZ we use the equation daz=AZ-AZ' (using approximation) and therefore AZ=AZ'+daz 0189 daz=azprime-azdoubleprime; 0190 del=elprime-eldoubleprime; 0191 0192 %Now we assume that the difference between the apparent and antenna 0193 %coordinates is small i.e that the same sort of correction would have been 0194 %found if we'd done the calculation on the apparent coordinates- defining 0195 %apparent coordinates as AZ and antenna coordinates as AZ' we have 0196 %daz = AZ-AZ' (not the difference and similarity (!!!) with the calculation 0197 %above 0198 % AZ=AZ'+daz 0199 0200 azapp = azprime+daz; 0201 elapp = elprime+del; 0202 0203 %So these are our first guesses for the apparent azimuth and elevation 0204 %values- they should be ok but lets just check by going from these values 0205 %back to the antenna values and comparing 0206 0207 [azprime2,elprime2]=calcAzEl2(azapp,elapp,d); 0208 0209 errazprime = (azprime-azprime2); 0210 errelprime = (elprime-elprime2); 0211 0212 maxazerr=max(abs(errazprime)); 0213 maxelerr=max(abs(errelprime)); 0214 vecazerr=[]; 0215 vecelerr=[]; 0216 a=isnan(daz); 0217 daz(a)=0; 0218 a=isnan(del); 0219 del(a)=0; 0220 %Now we iteratively adjust the Az (apparent coordinate) until it produces 0221 %the 'correct' AZ 0222 0223 alen=length(daz); 0224 0225 % for i=1:alen 0226 % if(~((isnan(azapp(i)) || isnan(elapp(i))))) 0227 % erraz(i)=1; 0228 % errel(i)=1; 0229 % dazt=azprime(i)-azdoubleprime(i); 0230 % delt=elprime(i)-eldoubleprime(i); 0231 % while(abs(erraz(i)) >0.00001 || abs(errel(i))>0.00001) 0232 % [azprime2(i),elprime2(i)]=calcAzEl3(azapp(i),elapp(i),d,i); 0233 % erraz(i) = azprime(i)-azprime2(i); 0234 % errel(i) = elprime(i) - elprime2(i); 0235 % daz(i)=dazt+erraz(i); 0236 % del(i)=delt+errel(i); 0237 % azapp(i) = wrapTo360(azprime(i)+daz(i)); 0238 % elapp(i) = wrapTo360(elprime(i)+del(i)); 0239 % end 0240 % end 0241 % end 0242 0243 while((maxazerr >0.00000001) || (maxelerr>0.00000001)) 0244 [azprime2,elprime2]=calcAzEl2(azapp,elapp,d); 0245 errazprime = (azprime - azprime2); 0246 errelprime = (elprime - elprime2); 0247 a=isnan(errazprime); 0248 errazprime(a)=0; 0249 a=isnan(errelprime); 0250 errelprime(a)=0; 0251 0252 maxazerr=max(abs(errazprime)); 0253 maxelerr=max(abs(errelprime)); 0254 daz = daz + errazprime; 0255 del = del + errelprime; 0256 azapp= wrapTo360(azprime+daz); 0257 elapp = wrapTo360(elprime+del); 0258 vecazerr=[vecazerr maxazerr]; 0259 vecelerr=[vecelerr maxelerr]; 0260 end 0261 0262 %might need some error checking here but otherwise 0263 0264 erraz=errazprime; 0265 errel=errelprime; 0266 aza=azapp; 0267 ela=elapp; 0268 0269 0270 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5 0271 function dmap = cleanRaDec_ovro_surf(dmap) 0272 0273 %this function takes a matrix in the form of dmap, and removes all the nan 0274 %rows basing this on columns 1 and 2 i.e assumes RA/DEC in columns 1 and 2 0275 %Used with the antenna to apparent coordinate conversion 0276 %CJC 0277 0278 a=isnan(dmap(:,1)); 0279 dmap(a,:)=[]; 0280 a=isnan(dmap(:,2)); 0281 dmap(a,:)=[];