% main.m
clc;
clear;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%timing
tic
%%%%%%%%%%%
%%%Continuously being updated
%%developed on P4 2.4Ghz laptop
%%750megs memory and Matlab 6.5 
%%%If it doesn't run on yours we may
%%%have to modify things
%%%send comments --> evr@u.washington.edu
%%%read Princeton camera *.SPE file
fid=fopen('630D13.SPE','r');
header=fread(fid,2050,'uint16'); %4100bytesHeader/2  
ImageMatrix=fread(fid,1024*1024,'uint16');
ImageMatrix=reshape(ImageMatrix,1024,1024);
fclose(fid);
ImageMatrix=double(ImageMatrix);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%plots%%%%%%%PLOT uncorrected *.SPE file******************
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%[X,Y]=meshgrid(1:1024,1:1024);
%figure(1);
%mesh(X,Y,ImageMatrix);
%view(90,90);% looks like Winview
%axis square;
%title('630D13.spe original');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%



%
%%%%%%%%%%%we want the center of the pattern%%%%%%%%%%%%%%%%
%there are about 8 (two are merged into one now and
%improvements can be made in locating these more accurately,etc.)
%defects fixed in
% the ccd camera .  These 8 hot spots are fixed in the 1024x1024
%grid but relative to the center of the pattern(xtry,ytry) these
%%will move.  
% we know that a circle is the shape of the pattern and we use that.
% assume we are at the center of the ImageMatrix xtry,ytry and send 
% IN rays every one degree searching for highest intensity store
%it and sort it dropping the lowest 
%Ndrop points of intensity.Find center and radius xc,yc and R.
%find the blocker arm relative to this center.
%search through the Ndrop points for the 
%longest continusous sequence
%%;that should be the angular coordinates 
%of the blocker arm relative
%to the assumed center xtry,ytry.  
%Construct the width of the arm and
%the ray through the center, 
%two lines distance width/2 parallel to 
%that ray and create mask of all
%points between those two lines
% fill in the blocker arm from points rotated about xc,yc.
%find the center Rx,Ry using least squares or weighted  1st moment
%find centroid of all points having r > R  make that new xtry,ytry
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%are there any corrections because (1,1) is ccd origin and not 0,0
%%no affects on circular symmetry because everything is relative
%%to xtry,ytry
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
RadiusBlocker=125;%assumed radius of blocker relative to xtry,ytry
xtry=512;%assumed center(to be iterated)
ytry=512;
%xold=512;% previous center(previous iteration)
%yold=512;
%%%%%%%%%%%%%%%%%%%No need to iterate
%%%%%%%%% by hand **No need to loop!!
%%%First and second iterations agree!!!!!!!!%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%begin iteration for finding center of pattern
%%%%%%%%
%%%%%%%%%%%%%%%%%location of hotspots%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%subject to almost certain change
%%%%%%%%%%%%
%nine hot spots  (x,y,z) of [lower left,upper right] 
%[(468,43,942),(514,6,984)] bright
%[(427,296,1416),(437,291,1128)] bright
% [(420,200,968),(438,196,524)] faint
% [(523,236,1672),(531,226,1562)] bright
% [(523,226, faint),(531,209,626)] faint combine with previous 
% [(475,654,959),(495,621,520)] bright
%[(461,716,1000),(465,712,1000)] faint guessing numbers
%[(397,898,1420),(427,878,1270)] bright
% [(574,907,785),(580,902,1015)] faint
MaskMatrix = ones(1024);

for x=468:514
    for y=6:43
        MaskMatrix(x,y)=0;
    end
end
for x=427:437
    for y=291:296
    MaskMatrix(x,y)=0;
    end
end
for x=420:438
    for y=196:200
        MaskMatrix(x,y)=0;
    end
end
for x=523:531
    for y=226:236
        MaskMatrix(x,y)=0;
    end
end
for x=523:531
   for y=209:226
    MaskMatrix(x,y)=0;
   end
end
for x=475:495
    for y=621:654
        MaskMatrix(x,y)=0;
    end
end
for x=461:465
    for y=712:716
       MaskMatrix(x,y)=0;
   end
end
for x=397:427
    for y=878:898
        MaskMatrix(x,y)=0;
    end
end
for x=574:580
    for y=902:907
        MaskMatrix(x,y)=0;
    end
end

%spy(MaskMatrix);
%view(90,90);





% TryMask  generated

%phi=angle to rotate MaskMatrix about (xtry,ytry)
phi = pi/2;
cosphi = cos(phi);
sinphi = sin(phi);
%TryMask will be zeros everywhere and ones where the hotspots are
TryMask = ones(1024) - MaskMatrix;
  % need to adjust for xtry,ytry
for x=1:1024
    for y=1:1024
        if(TryMask(x,y) == 1)
   xback = ceil(x*cosphi + y*sinphi + xtry*(1-cosphi)-ytry*sinphi);
   yback = ceil(-x*sinphi +y*cosphi +xtry*sinphi +ytry*(1-cosphi));
%check that xback and yback are within limits
if  xback >0 & xback <1025  
    if  yback >0 & yback <1025
        TryMask(x,y) = ImageMatrix(xback,yback);
    end%if yback
end% if xback    
end%if  TryMask  
    end%y
end%x
for x=2:1023
    for y=2:1023
        if(TryMask(x,y)==1)      %%%FIXed
TryMask(x,y)=ImageMatrix(2*xtry-x,2*ytry-y);
            %TryMask(x,y)=(TryMask(x+1,y)+TryMask(x-1,y)...
% +TryMask(x,y+1)+TryMask(x,y-1))/4;
        end%if
    end%y
end%x
        
%spy(TryMask);
%view(90,90);
TryMatrix = ImageMatrix .* MaskMatrix + TryMask;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%plot of 8 spots corrected on ccd%%%%%%%%%
%%%%%%%%%%%%
%figure(2);
%[X,Y]=meshgrid(1:1024,1:1024);
%mesh(X,Y,TryMatrix);
%axis square;
%title('automated correction of 8 hotspots on ccd');
%view(90,90);% looks like Winview
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%The blocker is off center so by useing a
%larger blocker radius we may help 
%RadiusBlocker may not be too critical because
% blocker very faint except at a few spots
%we find the max intensity and go outwards 
%from that --avoiding the blocker
%%%%%%%===%%%%%%%%%%%%%%%%===%%%%%%%%===%%%%%%%%%%%
%%%++++++%%%%+++++++++++%%%%%%%%%%%%%%%%+++%%%%%%%%
%360 x,y,z points store in xi,yi,zi
xi=zeros(1,360);
yi=zeros(1,360);
zi=zeros(1,360);

nrmax1 = min(xtry,1024-xtry) -7;%% could change 7 to something 
nrmax2 = min(ytry,1024-ytry) -7;% to avoid edge of ccd
nrmax=min(nrmax1,nrmax2);%nrmax should not run off ccd
nr=nrmax-RadiusBlocker;%number of radial points per angle
zith=zeros(1,nr);xith=zeros(1,nr);yith=zeros(1,nr);
%points of radius < RadiusBlocker and > nrmax are avoided

for ith=1:360
    for rth=nrmax:-1:RadiusBlocker % in toward center on a ray
        xith(rth)=ceil(xtry +rth*sin(ith*2*pi/360));
        %round up to nearest integer
        yith(rth)=ceil(ytry -rth*cos(ith*2*pi/360));
       
    zith(rth)=TryMatrix(xith(rth),yith(rth));     
    end%rth
    %find max zith and corresponding xith and yith and store
    %
  % if ( ith> 74 & ith< 89)% %%FOUND THE BLOCKER
 %   figure(ith)
 %     plot(zith,'*');
 % end
    [rv,mth]=sort(zith);%rv(nr)=maximum
    
    xi(ith)=xith(mth(nr));
    yi(ith)=yith(mth(nr));
   % xi(ith)=ceil(xtry +mth(nr)*sin(ith*2*pi/360));
    %yi(ith)=ceil(ytry -mth(nr)*cos(ith*2*pi/360));
    zi(ith)=rv(nr);
end% ith
%%%%% sort and drop lowest intensity
[zs,is]=sort(zi); %sorts intensities in ascending order
%zs=zi(is)
%assuming one maximum; sometime should fix it so
% the outermost maximum is used
Ndrop=160;%%number of points out of 360 to drop
%drop lowest Ndrop intensities 
%(this gets rid of straight part blocker) 
%and store Ncircle points 360-Ndrop
% xp,yp are Ncircle points on the 
%%circumference of max. intenity
%as viewed from approaching corners to center of ccd
Ncircle=360-Ndrop;

zp=zeros(1,Ncircle);
xp=zeros(1,Ncircle);
yp=zeros(1,Ncircle);


for ip=1:Ncircle  % Ncircle maximums
    zp(ip)=zs(is(ip+Ndrop));
    xp(ip)=xi(is(ip+Ndrop));
    yp(ip)=yi(is(ip+Ndrop));
end % ip
%figure(7);
%plot(xp,yp,'o') % plot the Ncircle  points on the circumference
%%%%%%%%%%%%%=======%%%%%%%%%%%%%%%%%%%%%%%%======%%%%%%%%%%%%%%%
%%%%%%%%%%%%%+++++++++++++%%%%%%%%%%%%%%%%%%+++%%%%%%%%%%%++%%%
%
%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%least squares%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find tentative center of these points using least squares
%and radius R
%   xp^2+yp^2 + a(1)*xp + a(2)*yp + a(3) = 0
nN=length(xp); xx=xp.*xp; yy=yp.*yp; xy=xp .* yp;
Als=[sum(xp) sum(yp) nN;sum(xy) sum(yy) sum(yp);sum(xx) sum(xy) sum(xp)];
Bls=[-sum(xx+yy); -sum(xx.*yp + yy .*yp);-sum(xx.*xp+xy.*yp)];
als=Als\Bls;
xc= -als(1)/2 ;  %center of circle
yc= -als(2)/2 ;
R=sqrt((als(1)^2 + als(2)^2)/4 -als(3)) ;   %radius of circle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% at this point we have rotated around the assumed center 
%to get 8 hotspots on ccd "smoothed"
%dropping the blocker stick from the circle of max. intensity
%using least squares on the Ncircle points that were not
%dropped finding (xc,yc) as center of radius R
% we will now extrapolate across the blocker
% arm .  Plot the intensities around the circle.
% 

%%%%%%%===%%%%%%%%%%%%%%%%===%%%%%%%%===%%%%%%%%%%%
%%%++++++%%%%+++++++++++%%%%%%%%%%%%%%%%+++%%%%%%%%
%360 x,y,z points store in xi,yi,zi
xi=zeros(1,360);
yi=zeros(1,360);
zi=zeros(1,360);


for ith=1:360
     rth=ceil(R) ;% around circle of radius R
        xi(ith)=ceil(xc +rth*sin(ith*2*pi/360));
        %round up to nearest integer
        yi(ith)=ceil(yc -rth*cos(ith*2*pi/360));
       
    zi(ith)=TryMatrix(xi(ith),yi(ith));     
      %[rv,mth]=sort(zi);%rv(360)=maximum
    
   % xi(ith)=xith(mth(nr));
   % yi(ith)=yith(mth(nr));
   % xi(ith)=ceil(xtry +mth(nr)*sin(ith*2*pi/360));
    %yi(ith)=ceil(ytry -mth(nr)*cos(ith*2*pi/360));
    %zi(ith)=rv(nr);
end% ith
%%%%% sort and drop lowest intensity
[zs,is]=sort(zi); %sorts intensities in ascending order
%zs=zi(is)
dy=diff(zs)./diff(is);
dx=[1:20];
dz=dy(1:20);ddz=diff(dy);
d2z=ddz(1:20);

%dy=[dy dy(359)];
%plot(dx,dz,dx,d2z)
[maxd2z,maxnum]=max(d2z);
[mind2z,minnum]=min(d2z);
%finding where second derivative changes sign
angularWidth=ceil(max(maxnum,minnum));
ais=is(1:angularWidth);
beginAngle=floor(min(ais)); % beginning angle of BLockerArm
   % on the circle of max. intensity
%endAngle =ceil(beginAngle + angularWidth -1 );%chenged this ***
if(angularWidth-floor((angularWidth/2))*2 )==1 %%angW is  ODD
    endAngle= beginAngle + angularWidth -1;
    midAngle= (beginAngle + endAngle)/2;
end% if odd angularWidth
if(angularWidth-floor((angularWidth/2))*2 )==0 %%angW is  even
    endAngle=beginAngle + angularWidth;
    midAngle=(beginAngle+endAngle)/2;
    
end% if even

  
%
%a point (x,y) relative to (xc,yc) on the circumference 
% of the circle of radius R is
%xR= xc +R*sin(theta*2*pi/360);
% y = yc -R*cos(theta*2*pi/360)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% made changes ****
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%do it wider %%%%%
%endAngle=2*endAngle -beginAngle +2;%double width +2 degrees
%beginAngle=beginAngle -2;% two degrees wider on top
RB=RadiusBlocker ; %%%%%%%%%%%%
%RB = RadiusBlocker;
MaskMatrix=ones(1024);
%MaskBlockArm=ones(1024);
%for theta=beginAngle:endAngle
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%start at blocker****************
%%%%%%%%%%%%%%%%%%%%%%
xB=ceil(xc +RB*sin((beginAngle)*2*pi/360));

%xR=ceil(xc +R*sin((beginAngle)*2*pi/360));
  %xR coord.on circle radius R
  midAngle=ceil((beginAngle + endAngle) /2);
 for x=xB:1024  
 %for x=round(xc):1024 %%%no effect%%
  for y=1:1024      
  %  two lines parallel to line of midAngle
 %if y1 <=y <=y2 mark it
   
   y1=(-x+xc+R*sin(beginAngle*2*pi/360))*cot(midAngle*2*pi/360);
  y1=ceil(y1+yc-R*cos(beginAngle*2*pi/360));% y coord. lies between
 
  y2=(-x+xc+R*sin(endAngle*2*pi/360))*cot(midAngle*2*pi/360);
  y2=ceil(y2+yc-R*cos(endAngle*2*pi/360));% y coord. lies between
 if (y1-40 <=y) & (y <= y2 + 40) % make wider %%%this WORKED!!!!!
 %if y1 <=y & y <= y2  % 
%MaskBlockArm(x,y)=0;
MaskMatrix(x,y)=0;
end % if
end %y
end %x
%spy(MaskMatrix);
%view(90,90);
%axis square;



%%%%%%%blocker%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%block%%%%%%%%
%%%%%%%blocker%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%block%%%%%%%%
%%%%%%%blocker%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%block%%%%%%%%
%TryMatrixB=TryMatrix;
%phiB=angle to rotate MaskBlockArm about (xc,yc)
phiB = - pi/4 ;
cosphi = cos(phiB);
sinphi = sin(phiB);
%TryMaskB will be zeros everywhere and ones where the hotspots are
TryMask=ones(1024) - MaskMatrix;
%TryMaskB = ones(1024) - MaskBlockArm;
  % need to adjust for xc,yc
for x=1:1024
    for y=1:1024
        if (TryMask(x,y) == 1)
        %if(TryMaskB(x,y) == 1)
   xback = ceil(x*cosphi + y*sinphi + xc*(1-cosphi)-yc*sinphi);
   yback = ceil(-x*sinphi +y*cosphi +xc*sinphi +yc*(1-cosphi));
%check that xback and yback are within limits
if  xback >0 & xback <1025 
    if  yback >0 & yback <1025
       TryMask(x,y) = TryMatrix(xback,yback);
        % TryMaskB(x,y) = TryMatrix(xback,yback);
    end%if yback
end% if xback    
end%if  TryMask  
    end%y
end%x
for x=2:1023
    for y=2:1023
       % if(TryMaskB(x,y)==1)
%TryMaskB(x,y)=(TryMaskB(x+1,y)+TryMaskB(x-1,y)+TryMaskB(x,y+1)
%%%%     +TryMaskB(x,y-1))/4;
  if(TryMask(x,y)==1)  %%%FIXed
   TryMask(x,y)=TryMatrix(2*xtry-x,2*ytry-y);
      
      %TryMask(x,y)=(TryMask(x+1,y)+TryMask(x-1,y)+...
   %TryMask(x,y+1)+TryMask(x,y-1))/4;
        end%if

        %end%if
    end%y
end%x
        
%spy(TryMask);
%view(90,90);
TryMatrix= TryMatrix .* MaskMatrix + TryMask;
%TryMatrixB = TryMatrix .* MaskBlockArm + TryMaskB;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%now  blocker arm outside of circle of maximum 
%intensity should be smoothed.  Let's see

%figure(16);
%[X,Y]=meshgrid(1:1024,1:1024);
%mesh(X,Y,TryMatrixB);
%axis square;
%view(90,90);% looks like Winview
%title('corrected hotspots and part of blocker arm');

%%%%%%%blocker%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%block%%%%%%%%
%%%%%%%blocker%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%block%%%%%%%%
%%%%%%%blocker%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%block%%%%%%%%

%%%%%%%%%%%%%%%%%Far Right%%%%%%%%%%%%%%%%
%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%FR%%%%%%%%%
MaskMatrix=ones(1024);
FR=985;
for y=1:1024
for x=FR:1024
    MaskMatrix(x,y)=0;
end%x
end%y
TryMask=ones(1024) - MaskMatrix;
%TryMaskB=ones(1024)-MaskMatrix;
phi=pi;% angle to rotate tryMaskB
cosphi = cos(phi);
sinphi = sin(phi);
%TryMaskB will be zeros everywhere and ones where the area of interest is

  % need to adjust for xc,yc
for x=1:1024
    for y=1:1024
        if(TryMask(x,y) == 1)
        %if(TryMaskB(x,y) == 1)
   xback = ceil(x*cosphi + y*sinphi + xc*(1-cosphi)-yc*sinphi);
   yback = ceil(-x*sinphi +y*cosphi +xc*sinphi +yc*(1-cosphi));
%check that xback and yback are within limits
if  xback >0 & xback <1025 
    if  yback >0 & yback <1025
        TryMask(x,y) = TryMatrix(xback,yback);
        %TryMaskB(x,y) = TryMatrix(xback,yback);
    end%if yback
end% if xback    
end%if  TryMask  
    end%y
end%x
for x=2:1023
    for y=2:1023
      if(TryMask(x,y)==1)
        %  if(TryMaskB(x,y)==1)
%TryMaskB(x,y)=(TryMaskB(x+1,y)+TryMaskB(x-1,y)+
%%       TryMaskB(x,y+1)+TryMaskB(x,y-1))/4;
TryMask(x,y)=TryMatrix(2*xtry-x,y);  %FIXed
%TryMask(x,y)=(TryMask(x+1,y)+TryMask(x-1,y)+TryMask(x,y+1)+TryMask(x,y-1))/4;
end%if
    end%y
end%x
        
%spy(TryMask);
%view(90,90);
TryMatrix=TryMatrix .* MaskMatrix + TryMask;
%TryMatrix = TryMatrixB .* MaskMatrix + TryMaskB;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%now  far right 
%intensity should be smoothed.  Let's see

%figure(17);
%[X,Y]=meshgrid(1:1024,1:1024);
%mesh(X,Y,TryMatrix);
%axis square;
%view(90,90);% looks like Winview
%title('630D13.spe corrected hotspots,part of blocker & far right');


%%%%%%%%%%%%%%%%%%%%%%%%%%%
% smootoh area inside blocker
%%




%%%%%%%%%%%%%%%%%Far Right%%%%%%%%%%%%%%%%
%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%FR%%%%%%%%%

%%%in retrospect I probably should have done this part
%%%%first and not last
%%%%%%%%%%%%%%%%%upper and lower edge%%%%%%%%%%%%%%%%
%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%upper and lower%%%%%%%%%
MaskMatrix=ones(1024);
%FR=985;
for x=1:1024
for y=1:12
    MaskMatrix(x,y)=0;
end%x
end%y
for x=1:1024
    for y=1020:1024
        MaskMatrix(x,y)=0;
    end
end

TryMask=ones(1024) - MaskMatrix;

phi=pi/2;% angle to rotate tryMask
cosphi = cos(phi);
sinphi = sin(phi);
%TryMaskB will be zeros everywhere and ones where the hotspots are

  % need to adjust for xc,yc
for x=1:1024
    for y=1:1024
        if(TryMask(x,y) == 1)
        
   xback = ceil(x*cosphi + y*sinphi + xc*(1-cosphi)-yc*sinphi);
   yback = ceil(-x*sinphi +y*cosphi +xc*sinphi +yc*(1-cosphi));
%check that xback and yback are within limits
if  xback >0 & xback <1025 
    if  yback >0 & yback <1025
        TryMask(x,y) = TryMatrix(xback,yback);
        
    end%if yback
end% if xback    
end%if  TryMask  
    end%y
end%x
for x=2:1023
    for y=2:1023
      if(TryMask(x,y)==1)  %%%%FIXed
 %if TryMask failed to grab the value what to use?
 % use a TryMask(x,y) with same radius
 TryMask(x,y)=TryMatrix(2*xtry-x,2*ytry-y);

%TryMask(x,y)=(TryMask(x+1,y)+TryMask(x-1,y)+TryMask(x,y+1)+TryMask(x,y-1))/4;
end%if
    end%y
end%x
        
%spy(TryMask);
%view(90,90);
TryMatrix=TryMatrix .* MaskMatrix + TryMask;







%%%%%%%%%%%%%%%%%%%%%%%%%
%%%begin centroid
%%%%%%%%%%%%%%
%%%%%
xcen=0;
ycen=0;
sumint=0;
for x=1:1024
    for y=1:1024
        %outside circle of maximum intensity
        if ((x-xc)*(x-xc) + (y-yc)*(y-yc)) >= R 
         intensity= TryMatrix(x,y);
            xcen = xcen + intensity * x;
         ycen = ycen + intensity * y;
         sumint=sumint + intensity;
     end% if >=R
 end %y
end%x
xcen=xcen/sumint;
ycen=ycen/sumint;
%center of the pattern%%%%%%%%%%%%%%%%%%%%%%center%%%%%%%%%%%%

xcen=xcen
ycen=ycen




%%%%%%%%%%%%%%%%%centroid found%%%%%%%%%%%%%%%%
%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%begin%%%%%%%%%%%%smooth center of blocker************
%%%%%%%%%%%%%%%%%%
%%  
%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%
MaskMatrix=ones(1024);%%al ones
TryMask=ones(1024)-MaskMatrix;%all zeros
minBlocker =100000000000;
for x=1:1024 %make smaller********************
for y=1:1024
    if (x-xcen)*(x-xcen) + (y-ycen)*(y-ycen) <= RB*RB
%find minimum value
if TryMatrix(x,y) < minBlocker
minBlocker=TryMatrix(x,y);
end % if

% TryMask(x,y)=(TryMask(x+1,y)+TryMask(x-1,y)+
%%%%   TryMask(x,y-1)+TryMask(x,y+1))/4;
end % if
end % y
end % x
RBB=RB + 30;% increase to 125 + 30
%%%%found minBlocker=min value of TryMatrix within the blocker
for x=1:1024 %make smaller********************
for y=1:1024
    if (x-xcen)*(x-xcen) + (y-ycen)*(y-ycen) <= RBB*RBB
      MaskMatrix(x,y)=0;
        TryMask(x,y)=minBlocker;
% TryMask(x,y)=(TryMask(x+1,y)+TryMask(x-1,y)+TryMask(x,y-1)+TryMask(x,y+1))/4;
end % if
end % y
end % x

TryMatrix=TryMatrix .* MaskMatrix + TryMask;

%%%%%%end%%%%%%smooth center of blocker %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%autoconvolution %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%
%nb=1024;
%B=TryMatrix;% TryMatrix 
%reverse conjugate
%B=conj(B(nb:-1:1,nb:-1:1));
%use power of 2 transform lengths
%nf=2^nextpow2(2*nb);
%at=fft2(B,nf,nf);
%bt=fft2(TryMatrix,nf,nf);
%multiply transform then inverse
%ZZ=ifft2(at.*bt);
%make real ouput for real input
%if ~any(any(imag(TryMatrix))) & ~any(any(imag(B)))
 %   ZZ = real(ZZ);
 %end
%trim
%ZZ(2*nb:nf,:)=[];
%ZZ(:,2*nb:nf)=[];





%[X,Y]=meshgrid(1:2047,1:2047);

%[X,Y]=meshgrid(1:1024,1:1024);
%figure(3);
%ZZ=conv2(Z,Z); this ran more than half hour & had to be killed
%mesh(X,Y,ZZ);
%view(90,90);
%axis square;
%figure(5);
%axis square;
%contourf(X,Y,ZZ);
%view(90,90);
%hold on; 
%shading flat;
%vf=(10^13)*linspace(1,3,50);
%[cf,hf]=contour(X,Y,ZZ,vf);
%clabel(cf,hf);
%colorbar;

%view(90,90);
%title('autoconvolution 630D13.SPE adjusted');
%print -djpeg99 Aout1.jpg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%now shift center of pattern to center of TryMatrix 
%and do autocorrelation using fourier transforms
cenMatrix=zeros(1024);

xcc=round(xcen-xtry);
ycc=round(ycen-ytry);
minOfTryMatrix=min(min(TryMatrix));
for x=1:1024
    for y=1:1024
 if (x-xcc)>0 & (x-xcc) < 1025 & (y-ycc) > 0 & (y-ycc) < 1025
                cenMatrix(x,y)=TryMatrix(x-xcc,y-ycc);
            else cenMatrix(x,y)=minOfTryMatrix;                
            end% if
    end%for y
end % for x

for x=2:1023
    for y=2:1023
      if cenMatrix(x,y) == 0
      

cenMatrix(x,y)=(cenMatrix(x+1,y)+cenMatrix(x-1,y)+cenMatrix(x,y+1));
cenMatrix(x,y)=(cenMatrix(x,y)+cenMatrix(x,y-1))/4  + minOfTryMatrix;
%%%%added minOfTryMatrix  was zero before
end%if
    end%y
end%x



%%%%%%%%%%%%%%%%%%%%%%%%finished centering
%%%%%%%%%%%%%%
%%%%%%%%Radial intensity%%%%%%%%%%%begin
%%%%%%
%This script is designed to calculate the Average Radial Intensity
%of a Data Set


%Define the center of the image determined from CenterofIntensity.m
%Ycenter = round(xcen);%517;%got x and y mixed up!!
%Xcenter = round(ycen);%503;% should use 512,512
%was off by 9
%fid=fopen('630D13.SPE','r');
%header=fread(fid,2050,'uint16'); %4100bytes/2
%ImMat=fread(fid,1024*1024,'uint16');
%Z=reshape(ImMat,1024,1024);
%fclose(fid);
%Z=double(Z);
%ImageMatrix=Z;
Ycenter = round(ytry); %512
Xcenter=round(xtry); %512
%Now lets find out how many bins there are:
Corner = zeros(1,4);
Corner(1) = sqrt((1-Xcenter)*(1-Xcenter)+(1-Ycenter)*(1-Ycenter));
Corner(2) = sqrt((1024-Xcenter)*(1024-Xcenter)+(1-Ycenter)*(1-Ycenter));
Corner(3) = sqrt((1-Xcenter)*(1-Xcenter)+(1024-Ycenter)*(1024-Ycenter));
Corner(4) = sqrt((1024-Xcenter)*(1024-Xcenter)+(1024-Ycenter)*(1024-Ycenter));
NumofBins = max(Corner);
NumofBins = fix(NumofBins); %truncate to integer


     
        
    %Now lets look at the radius of each pixel and place it in
    %its appropriate bin and then calculate the Average Intensity
    %for the radial bin.  To reduce the number of calculations
    %I will caluclate the radius for each pixel, round it down and then
    %declare it the Bin number. So each bin will be one pixel wide.
    
    BinIntensityVec = zeros(1,NumofBins);
    PixelsinBinVec = zeros(1,NumofBins);
    AvgBinIntensityVec = zeros(1,NumofBins);
    
    for m=1:1024
        for n=1:1024
            %Since Matlab doesn't like indicies=zero:
            if m==Xcenter&n==Ycenter
                continue
            end
            %Calculate Radius:
            Radius = sqrt((m-Xcenter)*(m-Xcenter)+(n-Ycenter)*(n-Ycenter));
         %For Bin: truncate to integer, ie Bin 4 runs from 4.000 to 4.999
            Bin=fix(Radius);
           % BinIntensityVec(Bin)=BinIntensityVec(Bin)+ImageMatrix(m,n);
            BinIntensityVec(Bin)=BinIntensityVec(Bin)+cenMatrix(m,n);
           
            PixelsinBinVec(Bin)=PixelsinBinVec(Bin)+1;
        end
    end
    
    %Now calculate the average Intensity in each Bin:
    FileArray = zeros(2,NumofBins);
    for j=1:NumofBins;
        AvgBinIntensityVec(j)=BinIntensityVec(j)/PixelsinBinVec(j);
        FileArray(1,j) = j;
        FileArray(2,j) = AvgBinIntensityVec(j);
    end
    
    %in terms of x and y
    x=FileArray(1,:);
    y=FileArray(2,:);
    
    %%%%%%%%%%%%%%%%%%%Levenberg-Marquardt non-linear least squares
    %%%%fitting to A*exp(-(x-B)^2)/2*C*C) + D*exp(-F*x)
    %%%%%%%%%%%%%FWHM=2.35 * C  ***p(3)=C D=p(4)  F=p(5)
    %%%%%%%%%%% B=p(2) A=p(1)/(p(3)*s2p)  s2p=sqrt(2*pi)
    %%%%%%%%%%%%A=max value B=center of gaussian 
    
    


% delete above and merge this into main
% function = gaussian plus an exponential (usually exponent
% is zero and the amplitude of exponential gives
%the base line  
% leasqr example/test
global verbose
%verbose(1)=0;%%%@@@@@@@@@@@I changed this
verbose(1) = 1; % normally 1 prints lots of intermediate results 
verbose (2) = 1;
%t=[200:1:700]';% assuming 1 pixel step
%t=[.1:.2:4.7]';
%t = [1:100]';%++++++x data
pin=[max(max(cenMatrix)),350,100,500,0]; %params 1 gauss 1 exponential
%pin=[.6 1 .3 3 .7]';% param guess 2 gaussians 5 param
t=(FileArray(1,:))';
data=(FileArray(2,:))';
%data=[1 6 28 62 130 165 117 81 36 24 42 29 ...
  %  41 45 52 40 25 30 23 10 3 6 2 2]';
%data=data./(.2*sum(data)); % normalize y data

%p = [1; 0.1]; 
%data = NLleasqrfunc (t, p);%++++++++++++
wt1 = (1 + 0 * t) ./ sqrt (data); % = 1 /sqrt of variances of data
%if (sscanf (version, "%f") < 4)
%%  disp ("Used 'rand'");
%  rand ("normal");
%++++++++
%  data = data + 0.05 * rand (100, 1) ./ wt1;   % 1 / wt1 = sqrt of var = standard deviation
%++++++++++ now supply y data

%else
 % disp ("Used 'randn'");
 % data = data + 0.05 * randn (100, 1) ./ wt1;
%endif;
options = [[0.01; 0.01;.01;.01;.01], [.8;.8;.8;.8;.8]];
%options = [[0.01; 0.01], [.8; .8]]; %+++++++
% Note by Thomas Walter :
%
% Using a step size of 1 to calculate the derivative is WRONG !!!!
% See numerical mathbooks why.
%[f1, p1, kvg1, iter1, corp1, covp1, covr1, stdresid1, Z1, r21] = ...
%  NLleasqr (t, data, [.8;.05], 'NLleasqrfunc', .001, 50, wt1, [1;1],  'NLleasqrdfdp', options);
% A derivative calculated from central differences needs:  step = 0.001...1.0e-8
% And a onesided needs:  step = 1.0e-5...1.0e-8 and may be still wrong

[f1, p1, kvg1, iter1, corp1, covp1, covr1, stdresid1, Z1, r21] = ...
  NLleasqr (t, data, pin, 'NLleasqrfunc', .001, 50, wt1, [0.001 ;.001;.001;.001;.001 ],  'NLleasqrdfdp', options);
% NLleasqr (t, data, [.8;.05], 'NLleasqrfunc', .001, 50, wt1, [0.001;0.001],  'NLdfdp', options);

% used in Xhw632
%[f1, p1, kvg1, iter1, corp1, covp1, covr1, stdresid1, Z1, r21] = ...
%  NLleasqr (t, data, pin, 'NLleasqrfunc', .001, 50, wt1, [0.001 ;.001;.001;.001;.001 ],  'NLdfdp', options);
%type Nlleasqr.m
%values of parameters returned are in vector p1 (p one)
height=(p1(1)/(sqrt(2*pi)*p1(3))) - p1(4) ;
width=2.3548*p1(3);
peakCenter=p1(2); 
concatT=strcat('cntr630D13.SPE Height= ',num2str(height),' Width= ',num2str(width),' peakCenter= ',num2str(peakCenter));
title(concatT);
legend('azimuthally averaged radial intensity','fitted to a gaussian',0);
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % figure(3);
  %  plot(x,y);% plot radial intensity
  %  plot(FileArray(1,:),FileArray(2,:))
    %title('radial intensity 630D13.SPE');%modify to include 
             %  the fit parameters
 % print -djpeg99 radialC630D13.jpg;%%%%%%%%%%@@@@@@@@@@@@@@@@
  
    %title('fitted to a1*exp(-(x-c1)/b1)^2 +a2*exp(-b2*x)');
 



%%%%%%%%end of Radial intensity%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%autoconvolution%%%%%%%%%%%%%%%%%%%%%%
%%%we really don't need it
%%% we need the fourier transform
% (or the inverse fourier transform ?)
nb=1024;
%nb=2048; %for Big
%%Big=zeros(2048);%##
%for x=1:1024
 %   for y=1:1024
  %      Big(512+x,512+y)=cenMatrix(x,y);
  % end
  %end
%B=Big;
%B=cenMatrix;% TryMatrix centered

%free memory

%free up memorey
clear MaskMatrix;
clear TryMatrix;
clear TryMask;
clear ImageMatrix;


%reverse conjugate
%B=conj(B(nb:-1:1,nb:-1:1));
%use power of 2 transform lengths
nf=2^nextpow2(2*nb);%##
%goodall ran out of memory sooner!!
%bt=fft2(Big,nf,nf);

%was able to get this far and ran out of memory
%reverse conjugate
%B=conj(B(nb:-1:1:nb:-1:1));

%at=fft2(B,nf,nf);

ZZcen=fft2(cenMatrix,nf,nf);%##
%bt=fft2(cenMatrix,nf,nf);
%multiply transform then inverse
%ZZcen=fft2(Big,nf,nf);

%ZZcen=ifft2(at.*bt);
%make real ouput for real input
%if ~any(any(imag(cenMatrix)))  %## & ~any(any(imag(B)))
 %  ZZcen = real(ZZcen);%##
 %end%##
%trim
%ZZcen(2*nb:nf,:)=[];%##
%ZZcen(:,2*nb:nf)=[];%##
%  clear cenMatrix;%##
%ZZcen=ifft(ZZcen,nf,nf);

%[X,Y]=meshgrid(1:4096,1:4096);
[X,Y]=meshgrid(1:2048,1:2048);%##

%[X,Y]=meshgrid(1:1024,1:1024);
%figure(3);
%ZZ=conv2(Z,Z); this ran more than half hour & had to be killed
%mesh(X,Y,ZZ);
%view(90,90);
%axis square;
%figure(6);
%axis square;
%contourf(X,Y,ZZ);
%view(90,90);
%hold on; 
%shading flat;
%vf=(10^13)*linspace(1,3,150);%changed 50 -> 150
%[cf,hf]=contour(X,Y,ZZcen,vf);

%clabel(cf,hf);
%colorbar;
%mesh(X,Y,ZZcen)
%view(90,90);

%title('autoconvolution 630D13.SPE adjusted and centered');
%figure(12);%##

%vf=(10^5)*linspace(1,3,50);%changed 150 -> 50
%meshc(X,Y,abs(ZZcen));

%ZZcen(1,1)=0;%this didn't help 
%ZZcen(2047,2047)=0;% still have spikes 
%ZZcen(2047,1)=0;%on four corners
%ZZcen(1,2047)=0;
ZZcen=fftshift(ZZcen);%##
ZZcen=ZZcen/(2047*2047);%##
ZZcen=abs(ZZcen);%## set ZZcen = to its absolute value


%for x=1000:1049  % take out centr hump
 %   for y=1000:1049% cnetered on 1024,1024
  %      ZZcen(x,y)=0;
  % end
  %end
  
RFT=100;%##
%RFT=120,60,100,55,50,45,40;35,30,25,20

  for x=1:2047%##
      for y=1:2047
          if ((x-1024)*(x-1024)+(y-1024)*(y-1024)) < RFT*RFT
              ZZcen(x,y)=0;
          end%if
      end%y
  end%x%##
  
%Thresh=5 .* 10^7.* zeros(2047);
%fThresh=find( ZZcen > Thresh); 
%ZZcen(fThresh) =0;% set those values of
%     ZZcen that are greater than Thresh
%     equalto zero and plot

mesh(X,Y,ZZcen);%##
view(90,90);%##
axis square;
%contour3(abs(ZZcen));
%waterfall(real(ZZcen));
%contour(X,Y,abs(ZZcen));
%[cf,hf]=contour3(X,Y,abs(ZZcen),vf);

%mx=max(max(ZZcen)); % computer went crazy

%mesh(X,Y,abs(ZZcen)./mx);

%vf=(10^9)*linspace(0.001,.2,50);
%[cf,hf]=contour(X,Y,abs(ZZcen),vf);% drove the computer carzy

%contour3(X,Y,ZZcen,vf);
%view(90,90);

%axis square;%##
%title('DC blocked modulus fourier Transform adjusted630D13.SPE/(2047*2047)'); 
% print -djpeg99 FT630D13.jpg

%title('autoconvolution Big630D13.SPE adjusted and centered');
%figure(13);
%mesh(X,Y,real(ZZcen));
%view(90,90);
%colorbar;
%axis square;
%title('real part of inverse fourier Transform Big630D13.SPE')
%print -djpeg99 Aout3.jpg


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

%now top and bottom 
%intensity should be smoothed.  Let's see
%%%%%%%%%%%final TryMatrix
%figure(18);
%[X,Y]=meshgrid(1:1024,1:1024);
%mesh(X,Y,TryMatrix);
%axis square;

%view(90,90);% looks like Winview

%title('630D13.spe corrected hotspots,part of blocker,right,top & bottom');
%print -djpeg99 Ain1.jpg
%%%%%%%%%%%%%%%end corrections hotspots,blocker,far right,top bottom%%%


%%%%%%%%%%%%%
%%%%I think everything after this is garbage%%%%%%
%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%to do%%%%%%%%%%%%%%
%assuming one maximum; sometime should fix it so
% the outermost maximum is used
%Ndrop=0;%%number of points out of 360 to drop
%Ncircle=360-Ndrop;

%zp=zeros(1,Ncircle);
%xp=zeros(1,Ncircle);
%yp=zeros(1,Ncircle); 
%for ip=1:Ncircle  % Ncircle maximums
%    zp(ip)=zs(is(ip+Ndrop));
%    xp(ip)=xi(is(ip+Ndrop));
%    yp(ip)=yi(is(ip+Ndrop));
%end % ip

%figure(7);
%plot(xp,yp,'o') % plot the Ncircle  points on the circumference
%view(90,90);
%axis square;
%title('around circle of radius R');
%%%%%%%%%%%%%=======%%%%%%%%%%%%%%%%%%%%%%%%======%%%%%%%%%%%%%%%
%the above could be streamlined
%%%%%%%%%%%%%+++++++++++++%%%%%%%%%%%%%%%%%%+++%%%%%%%%%%%++%%%

%sumint=0.;
%sumx=0.;
%sumy=0.;
%for ict=1:Ncircle % use the 200 angles
%    ith=is(ict);
%    sumrx=0.;
%    sumry=0.;
%   sumrint=0.;nb=1024;
%B=TryMatrix;% TryMatrix 
%reverse conjugate
%B=conj(B(nb:-1:1,nb:-1:1));
%use power of 2 transform lengths
%nf=2^nextpow2(2*nb);
%at=fft2(B,nf,nf);
%bt=fft2(TryMatrix,nf,nf);
%multiply transform then inverse
%ZZ=ifft2(at.*bt);
%make real ouput for real input
%if ~any(any(imag(TryMatrix))) & ~any(any(imag(B)))
 %   ZZ = real(ZZ);
 %end
%trim
%ZZ(2*nb:nf,:)=[];
%ZZ(:,2*nb:nf)=[];





%[X,Y]=meshgrid(1:2047,1:2047);

%[X,Y]=meshgrid(1:1024,1:1024);
%figure(3);
%ZZ=conv2(Z,Z); this ran more than half hour & had to be killed
%mesh(X,Y,ZZ);
%view(90,90);
%axis square;
%figure(5);
%axis square;
%contourf(X,Y,ZZ);
%view(90,90);
%hold on; 
%shading flat;
%vf=(10^13)*linspace(.75,1.75,20);
%[cf,hf]=contour(X,Y,ZZ,20);
%clabel(cf,hf);
%colorbar;

%view(90,90);
%title('autoconvolution 630D13.SPE adjusted');






%now top and bottom 
%intensity should be smoothed.  Let's see

%figure(18);
%[X,Y]=meshgrid(1:1024,1:1024);
%mesh(X,Y,TryMatrix);
%axis square;
%view(90,90);% looks like Winview
%title('630D13.spe corrected hotspots,part of blocker,right,top & bottom');

%%%%%%%%%%%%%%%end corrections hotspots,blocker,far right,top bottom%%%

%%%%%%%%%%%%%%%%%plotting of original%%%%%%%%%%%%%%%%%%%%%
%figure(19);
%[X,Y]=meshgrid(1:2047,1:2047);
%[X,Y]=meshgrid(1:1024,1:1024);
%mesh(X,Y,Big);
%mesh(X,Y,cenMatrix);
%axis square;
%view(90,90);% looks like Winview

%title( ' Big 630D13.SPE ');
%title('centered 630D13.spe corrected hotspots,blocker arm,right,top & bottom');
%print -djpeg99 Ain3.jpg % saves to file with highest quality jpg

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%to do%%%%%%%%%%%%%%
%assuming one maximum; sometime should fix it so
% the outermost maximum is used
%Ndrop=0;%%number of points out of 360 to drop
%Ncircle=360-Ndrop;

%zp=zeros(1,Ncircle);
%xp=zeros(1,Ncircle);
%yp=zeros(1,Ncircle); 
%for ip=1:Ncircle  % Ncircle maximums
%    zp(ip)=zs(is(ip+Ndrop));
%    xp(ip)=xi(is(ip+Ndrop));
%    yp(ip)=yi(is(ip+Ndrop));
%end % ip

%figure(7);
%plot(xp,yp,'o') % plot the Ncircle  points on the circumference
%view(90,90);
%axis square;
%title('around circle of radius R');
%%%%%%%%%%%%%=======%%%%%%%%%%%%%%%%%%%%%%%%======%%%%%%%%%%%%%%%
%the above could be streamlined
%%%%%%%%%%%%%+++++++++++++%%%%%%%%%%%%%%%%%%+++%%%%%%%%%%%++%%%

%sumint=0.;
%sumx=0.;
%sumy=0.;
%for ict=1:Ncircle % use the 200 angles

%   rthmax=ceil(sqrt(xp(ict)*xp(ict)+yp(ict)*yp(ict)));
%    for rth=nrmax:-1:rthmax % in toward center on a ray
       %%%%%%%%%%%%%%%%%%%%%%% begin here
 %       xith(rth)=ceil(xtry+rth*sin(ith*2*pi/360)); 
        %round up to nearest integer
 %       yith(rth)=ceil(ytry-rth*cos(ith*2*pi/360));
       
%    zith(rth)=ImageMatrix(xith(rth),yith(rth));     
%    sumrint=sumrint+rth*zith(rth);
    %replace this
%    sumrx = sumrx+zith(rth)*rth*rth*sin(ith*2*pi/360);
%    sumry = sumry-zith(rth)*rth*rth*cos(ith*2*pi/360);
%   end%rth
   
  
   
%  sumx=sumx+sumrx;
%   sumy=sumy+sumry;
%   sumint=sumint+sumrint;
% end% ith

%centerx=xtry + sumx/sumint;
%centery=ytry + sumy/sumint;


%it seems only natural to use square boxes instead of radius
% so call the previous centroidR.m and use x and y
% forget the theta=192-196 check this again  
%so 200 radial angles->less
% region of integration(summation) is in toward center on a ray until
% the maximum and omitting the blocker rays and ccd damaged rays
%the acts like it has a smaller area than 1024x1024
%seems to be rectangle   1<=x<=985  and 12<=y <=1020
%refine these numbers later
%       
 %%%%%%%%%%%%timing     

 
 
 toc