%This program will read the spike and EOD timing, together with swimming kinematics data from the output of
%the extractSpikesEODSwimDirSpeed.m, and will generate spatial plots of
%spikes both individuals and the spike rate plots
clear all
close all
exp_date=input('please type in the experiment date : ','s');
load(strcat(exp_date,'_Stage1_output.mat'));
 set(0,'defaultfigurecolor',[1 1 1])
%colormap jet
%gaussian smoothing factor
ks=0.1;
beta=-0.8;
fs=500;
%% Calculate the place where each spike ocurred. First interpolate the location of the fish at the time of each spike for each unit
for tn=1:length(hd_crct) %for all trials
    %% load the landmark locations data
    if filename{tn}(11)=='1' || filename{tn}(11)=='3' || filename{tn}(11)=='4' || filename{tn}(11)=='5'
        load(strcat('landmark_coordinates_',filename{tn}(1:8),'_t3.mat'));
        tankfilt=filter;
        clear filter
    else
        load(strcat('landmark_coordinates_',filename{tn},'.mat'));
        tankfilt=filter;
        clear filter
    end
    for b=1:length(xvalues{tn})
        pos_color(b,:)=[1-b/length(xvalues{tn}) 1-b/length(xvalues{tn}) 1-b/length(xvalues{tn})];
    end
    cnt=1;
    for i=1:length(nuraw_count{tn})
         for k=1:2
             Place_int_unit{tn,i,k}=interp1(xvalues{tn},hd_crct{tn}(:,k),nutimes{i,tn},'spline');
             if tn~=1
                Place_int_unit_tail{tn,i,k}=interp1(xvalues{tn},COM_crct{tn}(:,k),nutimes{i,tn},'spline');
             else
                Place_int_unit_tail{tn,i,k}=[];
             end

         end
         if nuraw_count{tn}(i)>10
            figure
             scatter(hd_crct{tn}(:,1),hd_crct{tn}(:,2),1,[0.8 0.8 0.8],'filled')
             hold on
             for b=1:nuraw_count{tn}(i)
                if Theta_interp{i,tn}(b)>95
                    mksize(b)=38;
                else
                    mksize(b)=8;
                end
                %find the closest time point in the head trajectory to the
                %time of the spike
                [m,n]=min(abs(xvalues{tn}-nutimes{i,tn}(b)));
                color_vect(b,:)=[1-pos_color(n,1) 0 0];
             end
             scatter(Place_int_unit{tn,i,1},Place_int_unit{tn,i,2},10,color_vect,'filled');
             %scatter(Place_int_unit_com{tn,i,1},Place_int_unit_com{tn,i,2},mksize',[0 1 0],'filled');
             clear color_vect
             clear mksize
             plotlm(tankfilt,tank_circle,hc,stp,ltp,sqp,sc,lc,filename{tn},1)  
             title(strcat(uname{i,tn},',DirInd =',num2str(DirPref{i,tn},'%4.2f\n')))
             axis off
             axis square
             cnt=cnt+1;
             %%saveas(gcf,strcat(filename{tn},'_',uname{i,tn},'_hd_SpikeSpatScat.pdf'))
            %%saveas(gcf,strcat(filename{tn},'_',uname{i,tn},'_hd_SpikeSpatScat.tif'))

         end
    end
    %do the same for the population spikes
    for k=1:2
        Place_int_pop{tn,k}=interp1(xvalues{tn},hd_crct{tn}(:,k),poptimes{tn},'spline');
    end
    clear color_vect
    for b=1:length(poptimes{tn})
        [m,n]=min(abs(xvalues{tn}-poptimes{tn}(b)));
        color_vect(b,:)=[1-pos_color(n,1) 0 0];
        if Theta_interp_POP{tn}(b)>95
           mksize(b)=38;
        else
          mksize(b)=8;
        end
    end
    h=figure;
    scatter(hd_crct{tn}(:,1),hd_crct{tn}(:,2),1,[0.8 0.8 0.8],'filled')
    hold on
    scatter(Place_int_pop{tn,1},Place_int_pop{tn,2},10,color_vect,'filled');
    plotlm(tankfilt,tank_circle,hc,stp,ltp,sqp,sc,lc,filename{tn},1)  
    title(strcat('Population',filename{tn},',DirInd =',num2str(DirPref_POP{tn},'%4.2f\n')))
    axis off
    axis square
    %%saveas(gcf,strcat(filename{tn},'hd_POP_SpikeSpatScat.pdf'))
    %%saveas(gcf,strcat(filename{tn},'hd_POP_SpikeSpatScat.tif'))
    
    %% calculate escans
    [vi, thresh2, vrAmpl, vnDur] = find_esac(EOD_crct{tn});
    %calculate the position of the fish at the time of escans
    X_int_escan=interp1(xvalues{tn},hd_crct{tn}(:,2),EOD_crct{tn}(vi),'spline');
    Y_int_escan=interp1(xvalues{tn},hd_crct{tn}(:,1),EOD_crct{tn}(vi),'spline');
    [vipa,vipl]=findpeaks((EOD_rate{tn}),'MinPeakHeight',mean(EOD_rate{tn}),'MinPeakWidth',5,'MinPeakProminence',0.5);

    figure
    plot(EOD_crct{tn},EOD_rate{tn})
    hold on
    plot(EOD_crct{tn}(vi),EOD_rate{tn}(vi),'*');
    plot(EOD_crct{tn}(vipl),EOD_rate{tn}(vipl),'ko')
    rate_thresh(tn)=prctile(EOD_rate{tn}(vi),97);
    RatePeak_index{tn}=vipl;
    Escan_index{tn}=vi;
    Escan_index_large{tn}=find(EOD_rate{tn}(vi)>rate_thresh(tn));
    Place_int_escan{tn}=interp1(xvalues{tn},hd_crct{tn},EOD_crct{tn}(vi));
    Place_int_rp{tn}=interp1(xvalues{tn},hd_crct{tn},EOD_crct{tn}(vipl));%locations of EOD rate peak

    Place_int_escan_large{tn}=interp1(xvalues{tn},hd_crct{tn},EOD_crct{tn}(vi(Escan_index_large{tn})));
    % find the largest peaks in the eod rate

    %%
    %scatter(Place_int_escan{tn}(:,1), Place_int_escan{tn}(:,2),5,[1 0 0],'filled');
    for b=1:length(Escan_index_large{tn})
            [m,n]=min(abs(xvalues{tn}-EOD_crct{tn}(vi(Escan_index_large{tn}(b)))));
            color_vect_escan(b,:)=[1-pos_color(n,1) 0 0];
    end
    figure(h.Number)
    scatter(Place_int_escan_large{tn}(:,1), Place_int_escan_large{tn}(:,2),10,color_vect_escan);

    axis square
    title('escan trajectory')
    clear pos_color
    clear color_vect
    clear color_vect_escan
    clear mksize
    clear vi
    %% calculate visit density. 
    %This is the number of times the fish visits a
    %given grid point. The visit density is insensitive to the variations in
    %swim speed since each grid passing was counted only once per visit (Jun et
    %al 2016)
    %use the location of the head interpolate to 100 Hz
    X_C{tn}=interp1(xvalues{tn},hd_crct{tn}(:,2),time{tn},'spline');
    Y_C{tn}=interp1(xvalues{tn},hd_crct{tn}(:,1),time{tn},'spline');
    y=0:80:1600;
    x=0:80:1200;
    resy=1600/80;
    resx=1200/80;
    visitd{tn}=NaN(resx,resy);
    visitt{tn}=zeros(resx,resy);
    timespent{tn}=zeros(resx,resy);

    temp=[0,0];
    p=1;
    %determin visit density for each bin
    for i=1:(resx)
        for j=1:(resy)
                temp1=((find(X_C{tn}>=x(i)&X_C{tn}<x(i+1)&Y_C{tn}>=y(j)&Y_C{tn}<y(j+1))));%use the 500 Hz interpolated fish location
                %do not include the times where the fish had stopped in the same location. Calculate only the first time it enters the location
                if length(temp1)==1
                    visitd{tn}(i,j)=1;
                elseif length(temp1)>1
                    visitd{tn}(i,j)=sum(diff(temp1)>1)+1;
                else
                    visitd{tn}(i,j)=NaN;
                end
                    timespent{tn}(i,j)=length(find(X_C{tn}>=x(i)&X_C{tn}<x(i+1)&Y_C{tn}>=y(j)&Y_C{tn}<y(j+1)))/fs;%time spent in each bin
                    clear temp2;
         end
    end
    

%% first calculate the total spike and escan count per bin 
    ndata=length(nuraw_count{tn});
    spkcnt{tn}=zeros(resx,resy,ndata);
    timespent_unit{tn}=zeros(resx,resy,ndata);
    occup{tn}=zeros(resx,resy,ndata);

    for w=1:ndata
        for k=1:length(Place_int_unit{tn,w,1})
            for i=1:resx
                for j=1:resy
                    if (Place_int_unit{tn,w,2}(k)>=x(i)& Place_int_unit{tn,w,2}(k)<x(i+1)& Place_int_unit{tn,w,1}(k)>=y(j)&Place_int_unit{tn,w,1}(k)<y(j+1))
                        spkcnt{tn}(i,j,w)=spkcnt{tn}(i,j,w)+1;
                    end
                end
            end
        end
        %delete isolated spike bins and call them NaN
        for ii=1:resx
            for jj=1:resy
                if spkcnt{tn}(ii,jj,w)>0
                    if ii>1 && jj>1 && ii<15
                        temp=nansum([spkcnt{tn}(ii-1,jj-1,w),spkcnt{tn}(ii,jj-1,w),spkcnt{tn}(ii+1,jj-1,w),spkcnt{tn}(ii-1,jj,w),spkcnt{tn}(ii+1,jj,w),spkcnt{tn}(ii-1,jj+1,w),spkcnt{tn}(ii,jj+1,w),spkcnt{tn}(ii+1,jj+1,w)]);
                    elseif ii==1 && jj>1
                        temp=nansum([spkcnt{tn}(ii,jj-1,w),spkcnt{tn}(ii+1,jj-1,w),spkcnt{tn}(ii+1,jj,w),spkcnt{tn}(ii,jj+1,w),spkcnt{tn}(ii+1,jj+1,w)]);
                    elseif ii>1 && jj==1
                        temp=nansum([spkcnt{tn}(ii-1,jj,w),spkcnt{tn}(ii+1,jj,w),spkcnt{tn}(ii-1,jj+1,w),spkcnt{tn}(ii,jj+1,w),spkcnt{tn}(ii+1,jj+1,w)]);
                    elseif ii==15
                        temp=nansum([spkcnt{tn}(ii-1,jj-1,w),spkcnt{tn}(ii,jj-1,w),spkcnt{tn}(ii-1,jj,w),spkcnt{tn}(ii-1,jj+1,w),spkcnt{tn}(ii,jj+1,w)]);
                    end
                    if temp<=0
                       spkcnt{tn}(ii,jj,w)=NaN;
                    end
                end
            end
        end
    end
    %if the visitdensity was less than 2, set timespent to NaN
    for ii=1:resx
        for jj=1:resy
            if visitd{tn}(ii,jj)<2
                visitd{tn}(ii,jj)=NaN;
                timespent{tn}(ii,jj)=NaN;
            end
        end
    end
    
    figure
    hold on
    pcolor([(visitd{tn}) NaN(resx,1);NaN(1,resy+1)])
    colormap(cool(20))
    brighten(beta)
    axis off
    shading flat
    title(strcat('visit density, ',num2str(max(max(visitd{tn})),'%4.2f\n')))
    plotlm(tankfilt,tank_circle,hc,stp,ltp,sqp,sc,lc,filename{tn},80)
    

%do the same thing for the population data
    spkcntpop{tn}=zeros(resx,resy);

    for k=1:length(Place_int_pop{tn,1})
            for i=1:resx
                for j=1:resy
                    if (Place_int_pop{tn,2}(k)>=x(i)&Place_int_pop{tn,2}(k)<x(i+1)& Place_int_pop{tn,1}(k)>=y(j)&Place_int_pop{tn,1}(k)<y(j+1))
                        spkcntpop{tn}(i,j)=spkcntpop{tn}(i,j)+1;
                    end
                end
            end
    end
 
    

    
    %% calculate visit density. Take only bins that contained more than 1 spike and more than 1 visit
    spkdens{tn}=NaN(resx,resy,ndata);
    spkdenspop{tn}=NaN(resx,resy);
    avefrbin{tn}=NaN(resx,resy,ndata);
    avefrbinpop{tn}=NaN(resx,resy);

    spkdens_smooth{tn}=NaN(resx,resy,ndata);
    spkdenspop_smooth{tn}=NaN(resx,resy);
    avefrbin_smooth{tn}=NaN(resx,resy,ndata);
    avefrbinpop_smooth{tn}=NaN(resx,resy);
    
    for w=1:ndata 
        %for each cell count only the timespent where spikecount is not NaN.
        timespent_unit{tn}(:,:,w)=timespent{tn}(:,:).*(~isnan(spkcnt{tn}(:,:,w)));
        total_time{tn}(w)=nansum(nansum(timespent_unit{tn}(:,:,w)));
        occup{tn}(:,:,w)=timespent_unit{tn}(:,:,w)/total_time{tn}(w);
        avefr{tn}(w)=nansum(nansum(spkcnt{tn}(:,:,w)))/total_time{tn}(w);
    
%     avefrbin_smooth{tn}(:,:,w)=imgaussfilt(spkcnt{tn}(:,:,w),ks,'FilterDomain','spatial')./imgaussfilt(timespent{tn},ks,'FilterDomain','spatial');
%     spkdens_smooth{tn}(:,:,w)=imgaussfilt(spkcnt{tn}(:,:,w),ks,'FilterDomain','spatial')./imgaussfilt(visitd{tn},ks,'FilterDomain','spatial');
    
        avefrbin{tn}(:,:,w)=spkcnt{tn}(:,:,w)./ timespent_unit{tn}(:,:,w);
        spkdens{tn}(:,:,w)=spkcnt{tn}(:,:,w)./visitd{tn};
    
    

   if nuraw_count{tn}(w)>10
        figure
        %imagesc(spkdens{tn}(:,:,w))
        
        pcolor([(avefrbin{tn}(:,:,w)) NaN(resx,1);NaN(1,resy+1)])
        colormap(cool(20))
        brighten(beta)
        axis off
        linfr=reshape((avefrbin{tn}(:,:,w)),[1,resx*resy]);
        %caxis([0,prctile(linfr,97)])
        shading flat
        %colorbar
        hold on
        title(strcat('spike frequency, ',uname{w,tn},', ',num2str(max(max(avefrbin{tn}(:,:,w))),'%4.2f\n')))
%         set(gca,'YDir','normal')
%         hold on
        [ii jj]=find(isnan(avefrbin{tn}(:,:,w))==1);
        [inside,on]=inpolygon(y(jj)+39,x(ii)+39,tank_circle(1,:)-39,tank_circle(2,:)-39);
        scatter(jj(inside)+0.5,ii(inside)+0.5, 's', 'MarkerEdgeColor','c','MarkerFaceColor','c')

        plotlm(tankfilt,tank_circle,hc,stp,ltp,sqp,sc,lc,filename{tn},80)
        %%%saveas(gcf,strcat(filename{tn},uname{w,tn},'Spikedens.pdf'))
        figure
        colormap(cool(20))
        brighten(beta)
        caxis([0,prctile(linfr,97)])
        colorbar
        %%%saveas(gcf,strcat(filename{tn},uname{w,tn},'Spikedens_cb.pdf'))

    end
end
%     avefrbinpop_smooth{tn}=imgaussfilt(spkcntpop{tn},ks,'FilterDomain','spatial')./imgaussfilt(timespent{tn},ks,'FilterDomain','spatial');
%     spkdenspop_smooth{tn}=imgaussfilt(spkcntpop{tn},ks,'FilterDomain','spatial')./imgaussfilt(visitd{tn},ks,'FilterDomain','spatial');
    avefrbinpop{tn}=spkcntpop{tn}./timespent{tn};
    spkdenspop{tn}=spkcntpop{tn}./visitd{tn};


    figure
    %imagesc(spkdenspop{tn}(:,:))
    pcolor([(avefrbinpop{tn}) NaN(resx,1);NaN(1,resy+1)])
        colormap(cool(20))
        brighten(beta)
    axis off

    shading flat
    %colorbar
    hold on
    title(strcat('population spike density, ',num2str(max(max(avefrbinpop{tn})),'%4.2f\n')))
    plotlm(tankfilt,tank_circle,hc,stp,ltp,sqp,sc,lc,filename{tn},80)
    %%saveas(gcf,strcat(filename{tn},uname{w,tn},'POPSpikedens.pdf'))

    %% calculate  escan density
    escancnt{tn}=zeros(resx,resy);
    escancntw{tn}=zeros(resx,resy);

    rpcnt{tn}=zeros(resx,resy);
    rpcntw{tn}=zeros(resx,resy);

for i=1:resx
    for j=1:resy
        for k=1:length(Place_int_escan{tn}(:,1))
            if (Place_int_escan{tn}(k,2)>=x(i)&Place_int_escan{tn}(k,2)<x(i+1)& Place_int_escan{tn}(k,1)>=y(j)&Place_int_escan{tn}(k,1)<y(j+1))
                        escancnt{tn}(i,j)=escancnt{tn}(i,j)+1;
                        escancntw{tn}(i,j)=escancnt{tn}(i,j)+EOD_rate{tn}(Escan_index{tn}(k));
            end
        end
        for k=1:length(Place_int_rp{tn}(:,1))
           if (Place_int_rp{tn}(k,2)>=x(i)&Place_int_rp{tn}(k,2)<x(i+1)& Place_int_rp{tn}(k,1)>=y(j)&Place_int_rp{tn}(k,1)<y(j+1))
                        rpcnt{tn}(i,j)=rpcnt{tn}(i,j)+1;
                        rpcntw{tn}(i,j)=rpcntw{tn}(i,j)+EOD_rate{tn}(RatePeak_index{tn}(k))
           end
        end

    end
end
escancnt_large{tn}=zeros(resx,resy);

for i=1:resx
    for j=1:resy
        for k=1:length(Place_int_escan_large{tn}(:,1))
        if (Place_int_escan_large{tn}(k,2)>=x(i)&Place_int_escan_large{tn}(k,2)<x(i+1)& Place_int_escan_large{tn}(k,1)>=y(j)&Place_int_escan_large{tn}(k,1)<y(j+1))
                    escancnt_large{tn}(i,j)=escancnt_large{tn}(i,j)+1;
        end
        end
    end
end

for i=1:resx
    for j=1:resy
        if isnan(visitd{tn}(i,j)) || visitd{tn}(i,j)<2
            escancnt_large{tn}(i,j)=NaN;
            escancnt{tn}(i,j)=NaN;
            escancntw{tn}(i,j)=NaN;
            rpcnt{tn}(i,j)=NaN;
            rpcntw{tn}(i,j)=NaN;
        end
    end
end

escandens{tn}=escancnt{tn}./visitd{tn};
escandensw{tn}=escancntw{tn}./visitd{tn};
escandens_large{tn}=escancnt_large{tn}./visitd{tn};
rpdens{tn}=rpcnt{tn}./visitd{tn};
rpdensw{tn}=rpcntw{tn}./visitd{tn};

escandens_smooth{tn}=imgaussfilt(escancnt{tn},ks,'FilterDomain','spatial')./imgaussfilt(visitd{tn},ks,'FilterDomain','spatial');
escandens_large_smooth{tn}=imgaussfilt(escancnt_large{tn},ks,'FilterDomain','spatial')./imgaussfilt(visitd{tn},ks,'FilterDomain','spatial');

%plot large escan density
figure
pcolor([(escandens_large{tn}) NaN(resx,1);NaN(1,resy+1)])
colormap(cool(20))
brighten(beta)
hold on
[ii jj]=find(isnan(escandens_large{tn})==1);
[inside,on]=inpolygon(y(jj)+39,x(ii)+39,tank_circle(1,:)-39,tank_circle(2,:)-39);
scatter(jj(inside)+0.5,ii(inside)+0.5, 's', 'MarkerEdgeColor','c','MarkerFaceColor','c')
plotlm(tankfilt,tank_circle,hc,stp,ltp,sqp,sc,lc,filename{tn},80)
axis off
shading flat
title(strcat('large escan density, ',num2str(max(max(escandens_large{tn})),'%4.2f\n')))
%%saveas(gcf,strcat(filename{tn},'LargeEscanDens.pdf'))
%plot all escans density
figure
pcolor([(escandens{tn}) NaN(resx,1);NaN(1,resy+1)])
colormap(cool(20))
brighten(beta)
hold on
[ii jj]=find(isnan(escandens{tn})==1);
[inside,on]=inpolygon(y(jj)+39,x(ii)+39,tank_circle(1,:)-39,tank_circle(2,:)-39);
scatter(jj(inside)+0.5,ii(inside)+0.5, 's', 'MarkerEdgeColor','c','MarkerFaceColor','c')
plotlm(tankfilt,tank_circle,hc,stp,ltp,sqp,sc,lc,filename{tn},80)
axis off
shading flat
title(strcat('all escan density, ',num2str(max(max(escandens_large{tn})),'%4.2f\n')))
%%saveas(gcf,strcat(filename{tn},'AllEscanDens.pdf'))

%plot weighted escan density

figure
pcolor([(escandensw{tn}) NaN(resx,1);NaN(1,resy+1)])
        colormap(cool(20))
        brighten(beta)
hold on
axis off
shading flat
title(strcat('weighted escan density, ',num2str(max(max(escandensw{tn})),'%4.2f\n')))
plotlm(tankfilt,tank_circle,hc,stp,ltp,sqp,sc,lc,filename{tn},80)

%plot peak rate density

figure
pcolor([(rpdens{tn}) NaN(resx,1);NaN(1,resy+1)])
        colormap(cool(20))
        brighten(beta)
hold on
axis off

%set(h, 'EdgeColor', 'b');
shading flat
title(strcat('rate peak density, ',num2str(max(max(rpdens{tn})),'%4.2f\n')))

% %colorbar
plotlm(tankfilt,tank_circle,hc,stp,ltp,sqp,sc,lc,filename{tn},80)

%%saveas(gcf,strcat(filename{tn},'_ratepeakdensity.pdf'))

%plot weighted peak rate density

figure
pcolor([(rpdensw{tn}) NaN(resx,1);NaN(1,resy+1)])
        colormap(cool(20))
        brighten(beta)
hold on
axis off

shading flat
title(strcat('weighted rate peak density, ',num2str(max(max(rpdensw{tn})),'%4.2f\n')))

% %colorbar
plotlm(tankfilt,tank_circle,hc,stp,ltp,sqp,sc,lc,filename{tn},80)
%%
for w=1:ndata
    %large escans
    temp_Escan_data=escandens_large{tn};
    temp_SPK_data=spkdens{tn}(:,:,w);
    Escan_data=temp_Escan_data(:);
    SPK_data=temp_SPK_data(:);
    
    [CC,p]=corrcoef(Escan_data(~isnan(Escan_data)& ~isnan(SPK_data)),SPK_data(~isnan(Escan_data)&~isnan(SPK_data)))
    SPK_Escan_corr{tn}(w,:)=[CC(1,2),p(1,2)];    
    %all escans
    temp_Escan_data=escandens{tn};
    AllEscan_data=temp_Escan_data(:);
    [CC,p]=corrcoef(AllEscan_data(~isnan(AllEscan_data)& ~isnan(SPK_data)),SPK_data(~isnan(AllEscan_data)&~isnan(SPK_data)))
    SPK_AllEscan_corr{tn}(w,:)=[CC(1,2),p(1,2)];    
    figure;plot(AllEscan_data,SPK_data,'*')
    %probability of spiking if there was no escans = probability of spiking
    %and no escans divided by probabilit of no escans
    NoEscan_SPK_Prob{tn}(w)=(sum(SPK_data(AllEscan_data==0 & ~isnan(SPK_data))>0)/sum(~isnan(SPK_data)))/(sum(AllEscan_data==0& ~isnan(SPK_data))>0)/sum(~isnan(AllEscan_data)&~isnan(SPK_data));
end
%%
close all
end
%%save(strcat(exp_date,'_Stage2_output.mat'),'NoEscan_SPK_Prob','SPK_AllEscan_corr','SPK_Escan_corr','EOD_crct','filename','time','nutimes','uname','RatePeak_index','Escan_index','Escan_index_large','poptimes','avefrbin','avefrbinpop','timespent','total_time','timespent_unit','occup','spkcnt','spkcntpop','avefr','visitd','escandens_large','escandensw','rpdensw','rate_thresh','Place_int_unit','Place_int_escan','Place_int_rp','Place_int_escan_large','Place_int_unit_tail');

%% Plot firing rate for each unit
%% Measure location specificity for each unit

function plotlm(tankfilt,tank_circle,hc,stp,ltp,sqp,sc,lc,filename,scale)
hold on
hold on
if scale ~=1
    scale1=scale;
    scale2=scale*1;
    a=1;
else
    scale1=scale;
    scale2=scale;
    a=0;
end

hold on
plot((tankfilt(:,1))/scale1+a,(1200-tankfilt(:,2))/scale2+a)
plot((tank_circle(1,:))/scale1+a,(1200-tank_circle(2,:))/scale2+a,'k')
if  filename(11)=='1'
        plot((hc(:,1))/scale1,(1200-hc(:,2))/scale2+a)

elseif filename(11)=='3'
        plot((hc(:,1))/scale1+a,(1200-hc(:,2))/scale2+a,':')
        plot((stp(:,1))/scale1+a,(1200-stp(:,2))/scale2+a);
        plot((ltp(:,1))/scale1+a,(1200-ltp(:,2))/scale2+a);
        plot((sqp(:,1))/scale1+a,(1200-sqp(:,2))/scale2+a);
        plot((sc(1,:))/scale1+a,(1200-sc(2,:))/scale2+a);
        plot((lc(1,:))/scale1+a,(1200-lc(2,:))/scale2+a);
elseif filename(11)=='4'
        plot(hc(:,1)/scale1+a,(1200-hc(:,2))/scale2+a,':')
        plot(ltp(:,1)/scale1+a,(1200-ltp(:,2))/scale2+a,':');

        plot(stp(:,1)/scale1+a,(1200-stp(:,2))/scale2+a);
        plot(sqp(:,1)/scale1+a,(1200-sqp(:,2))/scale2+a);
        plot(sc(1,:)/scale1+a,(1200-sc(2,:))/scale2+a);
        plot(lc(1,:)/scale1+a,(1200-lc(2,:))/scale2+a);
    
elseif filename(11)== '5'
        plot(hc(:,1)/scale1+a,(1200-hc(:,2))/scale2+a,':')
        plot(stp(:,1)/scale1+a,(1200-stp(:,2))/scale2+a);
        plot(ltp(:,1)/scale1+a,(1200-ltp(:,2))/scale2+a);
        plot(sqp(:,1)/scale1+a,(1200-sqp(:,2))/scale2+a);
        plot(sc(1,:)/scale1+a,(1200-sc(2,:))/scale2+a,':');
        plot(lc(1,:)/scale1+a,(1200-lc(2,:))/scale2+a);
elseif filename(11)== '6' || filename(11)== '7'
        plot(hc(:,1)/scale1+a,(1200-hc(:,2))/scale2+a,':')
        plot(stp(:,1)/scale1+a,(1200-stp(:,2))/scale2+a);
        plot(ltp(:,1)/scale1+a,(1200-ltp(:,2))/scale2+a);
        plot(sqp(:,1)/scale1+a,(1200-sqp(:,2))/scale2+a);
        plot(sc(1,:)/scale1+a,(1200-sc(2,:))/scale2+a,':');
        plot(lc(1,:)/scale1+a,(1200-lc(2,:))/scale2+a);
else
        plot(hc(:,1)/scale1+a,(1200-hc(:,2))/scale2+a,':')
        plot(stp(:,1)/scale1+a,(1200-stp(:,2))/scale2+a);
        plot(ltp(:,1)/scale1+a,(1200-ltp(:,2))/scale2+a);
        plot(sqp(:,1)/scale1+a,(1200-sqp(:,2))/scale2+a,':');
        plot(sc(1,:)/scale1+a,(1200-sc(2,:))/scale2+a,':');
        plot(lc(1,:)/scale1+a,(1200-lc(2,:))/scale2+a,':');
    
      
end
end

