function plotDSINC(N,p0,type)
% plot complex frequency response of 
% H(j*2*pi*N)=exp(j*(N-1)*pi*hz*p0).*(Sin(N*pi*hz*p0)./Sin(pi*hz*p0))
% for the Nyquist Band
% N is the number of sensors and 
% type is one of the following strings:
%    't'                      cycles per second,   p0=t0
%    'x'    x=1/lambda_0,     cycles per meter,    p0=x0
%    'phi'  phi=X/lambda_0,   wavelength response, p0=(x_0/lambda_0)*cos(phi)
%
% NOTES: Only the Nyquist Band is shown: -1/(2*p0) to 1/(2*p0)

warning off
figure

% set labels and plotting limits
if type=='t'
    label_x=['{\it{f}} [Hz]'];
    label_title=['Complex Freqency Response of {\it{H}}({\it{e^{j2\pi f t_0}}} )'];
    hz_limit=1/(2*p0);
    hz_step=.001/(2*p0);
    hz=-hz_limit:hz_step:hz_limit;    % create vector of cycles per unit (Hz)
else 
    if type == 'x'
        label_x=['1/{\it{\lambda}} [Hz]'];
        label_title=['Complex Freqency Response of {\it{H}}({\it{e^{j2\pi x_0 /\lambda}}} )'];
        hz_limit=1/(2*p0);
        hz_step=.001/(2*p0);
        hz=-hz_limit:hz_step:hz_limit;    % create vector of cycles per unit (Hz)    
    else 
        if type=='phi'
            label_x=['{\it{cos(\phi )}}'];
            label_title=['Complex Freqency Response of {\it{H}}({\it{e^{j2\pi N cos(\phi )}}} )'];
            hz_limit=1;
            hz_step=.001;
            hz=-hz_limit:hz_step:hz_limit;    % create vector of cycles per unit (Hz)    
        else
            error('Type must be either ''t'', ''x'', or ''phi''')
        end
    end
end


% Position Label
if type ~= 'phi'
    temp=1/(N*p0);
    point_label_x=temp-.55*temp;       
    point_label_y=-.03*N;
    point_label=num2str(temp);
else
    temp=1/(N*p0);
    point_label_x=temp-.1;       
    point_label_y=-.03*N;
    point_label=num2str(temp);
end

% Make Labels consist to six digits
temp2=length(point_label);
if temp2 < 6
    if round(temp)==temp             
        point_label=[point_label '.'];
        for i=temp2+2:6
            point_label=[point_label '0'];
        end
    else 
        for i=temp2+1:6
            point_label=[point_label '0'];
        end
    end
else
    if temp < 1000000
        point_label=point_label(1:6);
    end
end

% Complex Frequency Response
H_2pi_hz_p0=exp(-j*(N-1)*pi*hz*p0).*(sin(N*pi*hz*p0)./sin(pi*hz*p0));
                                  
% Calculate Magnitude of H(j*2*pi*f))                              
mag_H=abs(H_2pi_hz_p0);

% plot results
title(label_title)
hold on
axis square
axis([-hz_limit hz_limit -.1*N N])
plot(hz,mag_H)
plot([-hz_limit hz_limit],[0 0],'k')
xlabel(label_x), ylabel('Magnitude')
text(point_label_x,point_label_y,point_label)