function angplot(FINAL,HEIGHT,WIDTH,Xc,Yc)

% Puts intensity as a function of angle at defined radius.
% The result is a plot of the intensity along a circle
% centered on the crossing of the PSD-legs, at the radius
% the user is invited to input.
%
% The computation processes the average intensity on a window
% whose size can also be specified by the user.
%
% You may want to display the PSD before, to help you visualize
% and choose the appropriate radius.
%
% The main purpose of this program is to know the position of the PSD-legs.
% Once you know it, you can filter the PSD or run a radial scanning (radplot)
% along one of the legs.
%
% EXAMPLES:
%
% >> angplot(FINAL)
%    where FINAL is the PSD matrix.
%
% You can also specify the HEIGHT, the WIDTH, and the center coordinates Xc and Yc
% of your image. If you do so, type:
% >> angplot(FINAL,HEIGHT,WIDTH,Xc,Yc)
% Enter either 1 argument (i.e. FINAL), or 5.


if nargin==1
   [HEIGHT WIDTH]=size(FINAL);
   [Xc Yc]=center(FINAL);
end

if HEIGHT~=WIDTH
   disp(' WARNING: Your image isn''t square.')
   disp(' The result this program is going to return are WRONG...')
end

if nargin~=1 & nargin~=5
   error (' Number of expected input arguments: 1 or 5');
end


pi=3.14159265;

T=zeros(1,1024); % Intensity storage vector. Entries run through a circle at radius R0.
BUF=[]; % Temporary storage array

m=0;

while m~=3
   
   if m==0 | m==1
      iwid=input('\n Enter pixel window half-width for smoothing : ');
      sprintf(' The averaging window size is now %dx%d',2*iwid+1,2*iwid+1)
      sizeofwindow=(2*iwid+1)^2;
   end
   
   if m==0 | m==2
      R0=input(' Enter the radius (in pixels) of the circle \n you want to run through : ');
   end
   
   h=waitbar(0,'Angular scanning in progress...not that easy!');
   
   for I=1:1024
      waitbar(I/1024,h)
      TH=pi/512*(I-512); % Angle of study, from -pi to +pi
      CT=cos(TH);
      ST=sin(TH);
      X=Xc-R0*ST; % Real coordinates of study
      Y=Yc+R0*CT;
      M=fix(X); % Top left closest pixel from the real point of study
      N=fix(Y);
      
      % CAUTION: the origin is at (1,1), in the top left corner of the image.
      % That's why the expressions for X and Y are not those expected...
      % Thus, the image is described in the trigonometric way,
      % beginning on the right-hand side of the center.
      
      if 1<=M & M<=HEIGHT & 1<=N & N<=WIDTH % Stay away from edges
         for II=M:M+1 % Calculation around 4 pixels surrounding the real point
            for J=N:N+1
               LOCALSUM=0;
               for K=-iwid:iwid % Average of the window intensity
                  for L=-iwid:iwid
                     if 1<=II+K & II+K<=HEIGHT & 1<=J+L & J+L<=WIDTH % Prevents the window from overlapping out of the image
                        LOCALSUM=LOCALSUM+FINAL(II+K,J+L)/sizeofwindow;
                     end   
                  end
               end
               BUF(II,J)=LOCALSUM;
            end
         end
         
         % Now we'got 4 values : one on each pixel surrounding
         % the real point of study.
         % To synthetize those 4 values into 1, we will average them by 
         % weighing them with the area of the opposite little square
         % around the real point (the closer the pixel is from the point,
         % the more important it is in the average).
         
         P=X-M; % Lengths of the 4 squares
         Q=Y-N;
         PP=1-P;
         QQ=1-Q;
         T(I)=PP*QQ*BUF(M,N)+P*QQ*BUF(M+1,N)+Q*PP*BUF(M,N+1)+P*Q*BUF(M+1,N+1);
      end
   end
   
   close(h)
   
   figure;
   I=1:1024;
   ang=180/512*(I-512);
   plot(ang,T)
   axis on
   xlabel('Angle in degrees, 0 is on the right hand-side of the center')
   ylabel('Average intensity')
   title(sprintf('Angplot -- Intensity along the circle of radius %d pixels -- Averaging window size : %dx%d',R0,2*iwid+1,2*iwid+1))
   zoom on
   
   m=menu('WHAT WOULD YOU LIKE TO DO NOW?','Change window averaging size','Change the scanning radius','Quit or return to main menu');
   
end

return