apparentAzEl

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  function d = apparentAzEl(d, includeError)

  functions that retracts the pointing model and refraction corrections to
  go from encoder az/el to sky az/el.

  d - data structure
  
  creates new subfield:
    d.antenna0.servo.apparent, which is [Nby4], where each column
    corresponds to :
    [apparentAzimuth, apparentElevation, error in fit for Az, error in fit
    for El];


aza is the apparent azimuth
ela is the apparent elevation
erraz is the error between the antenna coordinate recorder in the data and
what we get when we use the apparent position to start with for refraction
and pointing calcualtions
similarly for erralt

  CJC

%%%%%%%%%%%%%%%%%%%%%%%%%%%
 ogk edit 1:
 Changed the condition in the while loop to
 Azimuth Error > threshold -> mod(Azimuth Error,360)>threshold
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  sjcm Aug 10, 2011:  modified inputs so we don't report back the error in
  az/el automatically (which i don't think actually gets used anywhere);
  memory problem.  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  Paddy Apr 11 2012: used wrap180 consistently throughout, avoids
  need for mod in threshold condition. Tidied up a bit.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function d = apparentAzEl(d, includeError)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 %  function d = apparentAzEl(d, includeError)
0006 %
0007 %  functions that retracts the pointing model and refraction corrections to
0008 %  go from encoder az/el to sky az/el.
0009 %
0010 %  d - data structure
0011 %
0012 %  creates new subfield:
0013 %    d.antenna0.servo.apparent, which is [Nby4], where each column
0014 %    corresponds to :
0015 %    [apparentAzimuth, apparentElevation, error in fit for Az, error in fit
0016 %    for El];
0017 %
0018 %
0019 %aza is the apparent azimuth
0020 %ela is the apparent elevation
0021 %erraz is the error between the antenna coordinate recorder in the data and
0022 %what we get when we use the apparent position to start with for refraction
0023 %and pointing calcualtions
0024 %similarly for erralt
0025 %
0026 %  CJC
0027 %
0028 %%%%%%%%%%%%%%%%%%%%%%%%%%%%
0029 % ogk edit 1:
0030 % Changed the condition in the while loop to
0031 % Azimuth Error > threshold -> mod(Azimuth Error,360)>threshold
0032 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0033 %  sjcm Aug 10, 2011:  modified inputs so we don't report back the error in
0034 %  az/el automatically (which i don't think actually gets used anywhere);
0035 %  memory problem.
0036 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0037 %  Paddy Apr 11 2012: used wrap180 consistently throughout, avoids
0038 %  need for mod in threshold condition. Tidied up a bit.
0039 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0040 
0041 if(nargin<2)
0042   includeError = 0;
0043 end
0044 
0045 %get the antenna coordinates that we start with
0046 azprime=d.antenna0.servo.az;
0047 elprime = d.antenna0.servo.el;
0048 maxitr=20;
0049 % if flag mask has been filled in, remove data with bits 0, 3 10 set
0050 % (corresponds to frame not received, tracker lacking status, and possible
0051 % servo timing issues).
0052 if(issubfield(d, 'flags', 'mask','fast'))
0053     bitmask = [0 3 10];
0054     indgood  = ~bitand(d.flags.mask.fast, sum(2.^bitmask)) & ...
0055     d.antenna0.servo.el>5 & d.antenna0.servo.el<88 & ...
0056     ~isnan(d.antenna0.servo.az) & ~isnan(d.antenna0.servo.el);
0057     dc = framecut(d, indgood);
0058 else
0059   indgood = true(size(azprime)) & d.antenna0.servo.el>5 & ...
0060       d.antenna0.servo.el<88 & ~isnan(d.antenna0.servo.el) & ~isnan(d.antenna0.servo.az);
0061   dc = framecut(d, indgood);
0062 end
0063 
0064 
0065 %lets just see what sort of correction we get at this antenna angle (AZ')-
0066 %if we assume the difference between AZ (apparent) and AZ'(antenna) is
0067 %small then we can use this number as a start approximation- here we have 3
0068 %coordinates Az (apparent AZ) Az' (antenna Az) AZ'' (meaningless azimuth
0069 %position calculated by applying a pointing and refraction correction to ...
0070 %    Az') and we
0071 %make the assumptio that Az-Az' is approximately Az'-Az''
0072 [azdoubleprime,eldoubleprime]=calcAzEl2(azprime(indgood),elprime(indgood),dc);
0073 
0074 %now we calculate what sort of offset we have, daz and del-note that this
0075 %offset is between AZ' and AZ''- ie daz=AZ'-AZ''- therefore to go back to
0076 %AZ we use the equation daz~=AZ-AZ' (using approximation) and therefore ...
0077 %    AZ=AZ'+daz
0078 
0079 daz=azprime(indgood)-azdoubleprime;
0080 del=elprime(indgood)-eldoubleprime;
0081 
0082 %Now we assume that the difference between the apparent and antenna
0083 %coordinates is small i.e that the same sort of correction would have been
0084 %found if we'd done the calculation on the apparent coordinates- defining
0085 %apparent coordinates as AZ and antenna coordinates as AZ' we have
0086 %daz = AZ-AZ' (not the difference and similarity (!!!) with the calculation
0087 %above
0088 % AZ=AZ'+daz
0089 
0090 azapp = wrap360(azprime(indgood)+daz);
0091 elapp = wrap360(elprime(indgood)+del);
0092 
0093 %So these are our first guesses for the apparent azimuth and elevation
0094 %values- they should be ok but lets just check by going from these values
0095 %back to the antenna values and comparing
0096 
0097 [azprime2,elprime2]=calcAzEl2(azapp,elapp,dc);
0098 
0099 errazprime = wrap180(azprime(indgood)-azprime2);
0100 errelprime = elprime(indgood)-elprime2;
0101 
0102 maxazerr=max(abs(errazprime));
0103 maxelerr=max(abs(errelprime));
0104 vecazerr=[];
0105 vecelerr=[];
0106 a=isnan(daz);
0107 daz(a)=0;
0108 a=isnan(del);
0109 del(a)=0;
0110 %Now we iteratively adjust the Az (apparent coordinate) until it produces
0111 %the 'correct' AZ
0112 
0113 %alen=length(daz);
0114 i=0;
0115 while(((maxazerr >0.0001) || (maxelerr>0.0001)) && (i<maxitr))
0116   i=i+1;
0117   disp(['Iteration (max): ', num2str(i),'(',num2str(maxitr),')']);  
0118   [azprime2,elprime2]=calcAzEl2(azapp,elapp,dc);
0119   errazprime = wrap180(azprime(indgood) - azprime2);
0120   errelprime = elprime(indgood) - elprime2;
0121   a=isnan(errazprime);
0122   errazprime(a)=0;
0123   a=isnan(errelprime);
0124   errelprime(a)=0;
0125   
0126   maxazerr=max(abs(errazprime));
0127   maxelerr=max(abs(errelprime));
0128   daz = daz + errazprime;
0129   del = del + errelprime;
0130   azapp = wrap360(azprime(indgood)+daz);
0131   elapp = wrap360(elprime(indgood)+del);
0132   vecazerr=[vecazerr maxazerr];
0133   vecelerr=[vecelerr maxelerr];
0134 %end
0135 end
0136 
0137 if i>=maxitr
0138     %Error handling in case the removal of the pointing correction isn't
0139     % working
0140     error('Apparent Az/El not converging')
0141 end
0142 %might need some error checking here but otherwise
0143 
0144 %erraz=errazprime;
0145 %errel=errelprime;
0146 %aza=azapp;
0147 %ela=elapp;
0148 totalIndices = 1:1:length(azprime);
0149 
0150 erraz = interp1(totalIndices(indgood), errazprime, totalIndices, ...
0151     'linear', 'extrap');
0152 errel = interp1(totalIndices(indgood), errelprime, totalIndices, ...
0153     'linear', 'extrap');
0154 aza   = interp1(totalIndices(indgood), azapp, totalIndices, ...
0155     'linear', 'extrap');
0156 ela   = interp1(totalIndices(indgood), elapp, totalIndices, ...
0157     'linear', 'extrap');
0158 
0159 
0160 
0161 
0162 
0163 
0164 if(includeError)
0165   d.antenna0.servo.apparent = [aza' ela' erraz' errel'];
0166 else
0167   d.antenna0.servo.apparent = [aza' ela']; 
0168 end
0169 
0170 return;