function ANGL=radplot(FINAL,HEIGHT,WIDTH,Xc,Yc)

% Puts intensity as a function of radius at defined angle.
% The result is a plot of the intensity along a leg of the PSD,
% centered on the crossing of the PSD-legs, at the angle
% the user is invited to input (this angle is set as an output,
% in case it is needed afterwards).
%
% The computation processes the average intensity on a window
% whose size can also be specified by the user.
%
% You may want to run the program angplot before, to know the accurate
% position of the PSD-legs.
%
% EXAMPLES:
%
% >> radplot(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.
%
% To get the angle of the scanned leg and store it in THETA, type:
% >> THETA=radplot(...);

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

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


pi=3.14159265;

T=zeros(1,2049); % Intensity storage vector. Entries run along a line at angle TH0.
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 window averaging size is now %dx%d',2*iwid+1,2*iwid+1)
      sizeofwindow=(2*iwid+1)^2;
   end
   
   if m==0 | m==2
      TH0=input(' Enter the angle (in degrees) of the leg \n you want to run along : ');
      ANGL=TH0; % Just for the convenience. ANGL for the output, TH0 in the program.
   end
   
   h=waitbar(0,'Radial scanning in progress...quite tough!');
   
   for I=-1024:1024
      waitbar((I+1024)/2048,h)
      R=sqrt((HEIGHT/2)^2+(WIDTH/2)^2)/1024*I; % Radius of study
      CT=cos(TH0*2*pi/360);
      ST=sin(TH0*2*pi/360);
      X=Xc-R*ST; % Real coordinates of study
      Y=Yc+R*CT;
      M=fix(X); % Top left closest pixel from the real point of study
      N=fix(Y);
      
      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+1024)=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=-1024:1024; rad=(sqrt((HEIGHT/2)^2+(WIDTH/2)^2)/1024)*I;
   plot(rad,T)
   axis on
   xlabel('Radius in pixels from the center')
   ylabel('Average intensity')
   title(sprintf('Radplot -- Intensity along the line of angle %d degrees -- Averaging window size : %dx%d',TH0,2*iwid+1,2*iwid+1))
   zoom on
   
   m=menu('WHAT WOULD YOU LIKE TO DO NOW?','Change window averaging size','Change the scanning angle','Quit or return to main menu');
   
end

return