temdiff.m

Calculate the integrated intensity of electron diffraction patterns.

function Distribution
% This M file is used to calculate the integrated intensity of Electron
% Diffraction patterns, I(r), where r is the radius of diffraction rings. 3/7/2006
[filename, pathname] = uigetfile('*.*', 'Read data file');
yc=input('The X value of center spot   ');
xc=input('The Y value of center spot   ');
rmin=input('The minimum radius   ');
rmax=input('The maximum radius   ');
dr=input('The step of radius (Must be an interger, minimum is 1 pixel!)    ');
dr=dr*2;
dtheta=input('The step of polar angle along the diffraction rings   ');
pi=3.1415926;
h=4.1356692e-15;
c0=2.99792458e8;
m0=0.51099906/c0/c0*1e6;
e0=1;
scale=input('The scale of image: The length in mm per pixel; If you do not know the scale, input 0---Hitachi; 1---Techai  ');
if scale==0
else
    mode=input('Type of microscopy: 1---Hitachi; 2---Tecnai ');
    if mode==1
    len=input('The length of camera in mm   ');
    else
        scale1=input('The scale of diffraction Ring (nm-1)   ');
        scale2=input('The scale bar length in pixel:  ');
        scale3=scale2/scale1;
    end
end
U0=input('The accelerated voltage  kV   ');
U0=U0*1000;
lamda=h/sqrt(2*m0*e0*U0*(1+e0*U0/2/m0/c0/c0))*1e9;
[filename1, pathname] = uiputfile('*.dat', 'Save data file');
data=dlmread(filename);
data=data';
[MaxX, MaxY]=size(data);
for Ir=rmin:dr:rmax
    ResultsN(Ir,1)=Ir;
    ResultsN(Ir,2)=0;
end
for IIr=rmin:dr:rmax
    r1=IIr-dr/2;
    r2=IIr+dr/2;
    for I=1:360/dtheta+1
        result(I,1)=(I-1)*dtheta;
        result(I,2)=0;
    end
    for II=xc-r2:xc
       if II>0
            for JJ=yc-r2:yc+r2
                if JJ>0
                    if JJ<MaxX
                        x=II+1;
                        y=JJ+1;
                        r=sqrt((x-xc)^2+(y-yc)^2);
                        if r>=r1
                           if r<=r2
                               theta=acos((y-yc)/r)*180/pi;
                               N=round(theta/dtheta)+1;
                               result(N,2)=result(N,2)+data(y,x);
                           end
                        end
                    end
                end
            end
       end
    end
    for III=xc+1:xc+r2
        if III<MaxY
            for JJJ=yc-r2:yc+r2
                if JJJ>0
                    if JJJ<MaxX
                        x=III+1;
                        y=JJJ+1;
                        r=sqrt((x-xc)^2+(y-yc)^2);
                        if r>=r1
                           if r<=r2
                               theta=acos((-y+yc)/r)*180/pi+180;
                               N=round(theta/dtheta)+1;
                               result(N,2)=result(N,2)+data(y,x);
                           end
                        end
                    end
                end
            end
        end
    end
    result(1,2)=result(1,2)+result(360/dtheta+1,2);
    result(360/dtheta+1,2)=result(1,2);
    for I2=1:360/dtheta+1
        ResultsN(IIr,2)=ResultsN(IIr,2)+result(I2,2);
    end
end
f2=fopen(filename1,'w');
if scale==0
    for I3=rmin:dr:rmax
     fprintf(f2,'%d  %f \n', ResultsN(I3,1),ResultsN(I3,2));
    end
    I4=rmin:dr:rmax;
    plot(ResultsN(I4,1), ResultsN(I4,2), 'r-');
    title('Intensity vs d-spacing');
    xlabel('d (pixel)');
    ylabel('Intensity');
else
    if mode==1
    for I3=rmin:dr:rmax
     fprintf(f2,'%d  %f \n', lamda*len/scale/ResultsN(I3,1),ResultsN(I3,2));
     res(I3,1)=lamda*len/scale/ResultsN(I3,1);
     res(I3,2)=ResultsN(I3,2);
    end
    else
        for I3=rmin:dr:rmax
        fprintf(f2,'%d  %f \n', scale3/ResultsN(I3,1),ResultsN(I3,2));
        res(I3,1)=scale3/ResultsN(I3,1);
        res(I3,2)=ResultsN(I3,2);
        end
    end
    I4=rmin:dr:rmax;
    plot(res(I4,1),res(I4,2),'r-');
    title('Intensity vs d-spacing');
    xlabel('d (nm)');
    ylabel('Intensity');         
end
status=fclose(f2);