%Script that generates Figure 8C from Hasselmo, M.E. (2009) A model of
%episodic memory: Mental time travel along encoded trajectories using grid
%cells.  Neurobiol. Learning and Memory 92: 559-573.

%Script written by
%Michael E. Hasselmo, Boston University, copyright 2007, 2008.
 
close all
clear all

tic

%TIME
datainput=8; %datainput=0=random, 1 - Hafting, datainput=4  = circular, 5=circular with stopping, 6=mapped route 
%7=shorter map with stim, 8 overlapping with stim

fullT=300; %START.    
kplot=60; k=0; oldt=1; %Plotting
kplace=1; kp=kplace; %Interval of place cell formation. 

%HEAD DIRECTION CELL PARAMETERS
dt=0.02;
h=[0; pi/3; 2*pi/3; pi; 4*pi/3; 5*pi/3]+pi/5;
HHst=[cos(h) sin(h)];

hd=zeros(length(h),fullT); %Head direction cell activity at each time.
HDnum=length(h);
arcnum=1;
numP=1; %Total=numP^2 (e.g. numP=10 gives 100 place cells)

%RUN REPEATED TESTS WITH DIFFERENT PARAMETERS
numPsizesTot=1; %Number of different sizes of numP or other variable. Get memory errors when numP>36.
numPtestsTot=1; % 3%Number of individual tests at a single variable value.
numCircRun=zeros(numPsizesTot,numPtestsTot+1);
numPlaceToPlot=1; %This sets number of place cells that will be plotted in place cell plots.
exdistR=zeros(numPsizesTot,numPtestsTot+1);
stopWhile=0;
strengHD=1.0;


%FIGURE FOR PLOTTING TRAJECTORIES
figSize=numPsizesTot;
figure('name','x y enc','Numbertitle','off','Position',[360 700 300 300]); ax14b=gca;
figure('name','x y ret','Numbertitle','off','Position',[1010 700 300 300]); ax14=gca;
axes(ax14);

rc=[1 0 0; 1 1 0; 1 0 1; 1 1 1; 0 1 1; 0 1 0.5; 0 1 0; 1 1 0; 1 0.5 0; 1 0 0];
colorSpikes=0;
figure('name','x,y,speed vs time','Numbertitle','off','Position',[50 700 300 100+70]); ax111=gca;
figure('name','x,y,speed vs time retrieval','Numbertitle','off','Position',[700 700 300 100+70]); ax112=gca;
figure('name','HD angle vs time','Numbertitle','off','Position',[360 400 300 100+70]); ax113=gca;
figure('name','HD angle vs time retrieval','Numbertitle','off','Position',[1010 400 300 100+70]); ax114=gca;
figure('name','Place Enc v time','Numbertitle','off','Position',[50 400 300 100+70*figSize]); ax16b=gca;
figure('name','Place Ret v time','Numbertitle','off','Position',[700 400 300 100+70*figSize]); ax21b=gca;
figure('name','HD encoding vs time','Numbertitle','off','Position',[50 100 300 100+70]); ax18b=gca;
figure('name','HD retrieval vs time','Numbertitle','off','Position',[700 100 300 100+70*figSize]); ax19b=gca;
figure('name','HD vectors Enc','Numbertitle','off','Position',[360 100 160 100+70]); ax18=gca;
figure('name','HD vectors Ret','Numbertitle','off','Position',[1010 100 160 100+70*figSize]); ax19=gca; 

figure('name','dist ret-enc','Numbertitle','off','Position',[400 50 100*numPsizesTot 50*numPtestsTot]); ax22x=gca;

for numPsizes=1:numPsizesTot;
    numP=1;  %Sets the number of place cells.  Get retrieval of nonoverlapping components that jumps ahead to end of trajectory.
    arcnum=50;   %Sets the number of arc length cells.  Get good retrieval of overlapping trajectory.

    if numP==0; numP=1; end
    
    for numPtests=1:numPtestsTot;
            numCircRun(numPsizes,numPtests+1)=1;
            placeSpike2=zeros(numP^2,fullT);
            placeSpike3=zeros(numP^2,fullT);
            plotrphi=zeros(length(h)/2,fullT);
            hd=zeros(length(h),fullT);
            hd3=zeros(length(h),fullT);
            fpsR=[];
            fpsE=[];
            fpsP=[];
            fpsPR=[];
            fpsHDeee=[];
            fpsHDrrr=[];
            spikehd=[];
            spikehd3=[];



%CREATE TRAJECTORIES - 
speed=5; revs=0.5; mom=0.99; %Rat movement
xv=zeros(2,fullT); 
dxv=zeros(2,fullT); %Position vectors.

minT=1;

if datainput==1;
    minT=8; %1208 
load ..\Hafting_Fig2c_Trial2.mat
dt=0.02; sumspeed=0;
xv(1,1:fullT+1)=pos_x(minT:minT+fullT,1); %Take position from data file pos_x at step t+minT.
xv(2,1:fullT+1)=pos_y(minT:minT+fullT,1); %Take position from data file pos_y at step t
dxv(1,1:fullT+1)=(pos_x(minT:minT+fullT,1)-pos_x(minT-1:minT+fullT-1,1))/dt;
dxv(2,1:fullT+1)=(pos_y(minT:minT+fullT,1)-pos_y(minT-1:minT+fullT-1,1))/dt; %Velocity vector over time.
        for t=1:fullT
            espeed=sqrt(dxv(1,t+1)^2+dxv(2,t+1)^2);
            sumspeed=sumspeed+speed;
        end
        sumspeed*dt %Print total distance covered.
        fullT*dt     %Print total time run.
        avespeed=sumspeed/fullT %Print average speed.
        xmax=100; xmin=-100;  ymax=100; ymin=-100;

elseif datainput==6; %MAP OF MOVEMENT
        dt=0.02; sumspeed=0;
    numCircRun=[floor(strengHD*4)-1 floor(strengHD*4)-1]; %Stops cycling through different runs for "while" command
    minT=1; %For random trajectory, must set minT=2 so that t-1 is at index 1.
    xmax=100; xmin=-100; ymax=180; ymin=-20; %Environment size.
    xv(:,1)=[30,30]; 
    tdur=0; t=1;
    tNdur=0;
    Map=[300   0  10; 
         300 -10   0;
         400   0  10;
         200  10   0;
         200   0 -10;
         500 -10   0];
    for mm=1:length(Map(:,1));
    dxv(:,t)=Map(mm, 2:3)';
    tdur=tNdur; tNdur=Map(mm,1)+tdur;
    for t=tdur+1:tNdur+1; dxv(:,t+1)=dxv(:,t); xv(:,t+1) = xv(:,t) + dt*dxv(:,t+1); end
    end
    pos_x=xv(1,:);
    pos_y=xv(2,:);
    fullT=length(pos_y)-1;
    espeed=zeros(1,fullT+1);
        for t=1:fullT
           espeed(1,t)=sqrt(dxv(1,t+1)^2+dxv(2,t+1)^2);
        end
elseif datainput==7; %MAP OF MOVEMENT
        dt=0.02; sumspeed=0;
    numCircRun=[floor(strengHD*4)-1 floor(strengHD*4)-1]; %Stops cycling through different runs for "while" command
    minT=1; %For random trajectory, must set minT=2 so that t-1 is at index 1.
    xmax=100; xmin=-100; ymax=100; ymin=-100; %Environment size.
    xv(:,1)=[30,-80]; 
    tdur=0; t=1;
    tNdur=0;
     Map=[300   0  10; 
         300 -10   0;
         400   0  10;
         400  10   0];
     Stim=[40 -30 1;
           -40 -20 2;
           -40 40 3]
     StimVect=[1 0 0; %One column for each STIMULUS!
               0 1 0;
               0 0 1]
    for mm=1:length(Map(:,1));
    dxv(:,t)=Map(mm, 2:3)';
    tdur=tNdur; tNdur=Map(mm,1)+tdur;
    for t=tdur+1:tNdur+1; dxv(:,t+1)=dxv(:,t); xv(:,t+1) = xv(:,t) + dt*dxv(:,t+1); end
    end
    pos_x=xv(1,:);
    pos_y=xv(2,:);
    fullT=length(pos_y)-1;
    espeed=zeros(1,fullT+1);
        for t=1:fullT
           espeed(1,t)=sqrt(dxv(1,t+1)^2+dxv(2,t+1)^2);
        end  
elseif datainput==8; %OVERLAPPING MAP OF MOVEMENT
        dt=0.02; sumspeed=0;
    numCircRun=[floor(strengHD*4)-1 floor(strengHD*4)-1]; %Stops cycling through different runs for "while" command
    minT=1; %For random trajectory, must set minT=2 so that t-1 is at index 1.
    xmax=100; xmin=-100; ymax=100; ymin=-100; %Environment size.
    xv(:,1)=[30,-80]; 
    tdur=0; t=1;
    tNdur=0;
        Map=[300   0  10; 
         300 -10   0;
         400   0  10;
         200  10   0;
         200   0 -10;
         500 -10   0];
     Stim=[40 -30 1;
           -40 -20 2;
           -40 40 3;
           20 40 4;
           -80 10 5]
     StimVect=[1 0 0 0 0; %One column for each STIMULUS!
               0 1 0 0 0;
               0 0 1 0 0;
               0 0 0 1 0;
               0 0 0 0 1]
    for mm=1:length(Map(:,1));
    dxv(:,t)=Map(mm, 2:3)';
    tdur=tNdur; tNdur=Map(mm,1)+tdur;
    for t=tdur+1:tNdur+1; dxv(:,t+1)=dxv(:,t); xv(:,t+1) = xv(:,t) + dt*dxv(:,t+1); end
    end
    pos_x=xv(1,:);
    pos_y=xv(2,:);
    fullT=length(pos_y)-1;
    espeed=zeros(1,fullT+1);
        for t=1:fullT
           espeed(1,t)=sqrt(dxv(1,t+1)^2+dxv(2,t+1)^2);
        end  
elseif datainput==0; %RANDOM MOVEMENT
    dt=0.02; sumspeed=0;
    numCircRun=[floor(strengHD*4)-1 floor(strengHD*4)-1]; %Stops cycling through different runs for "while" command
    minT=1; %For random trajectory, must set minT=2 so that t-1 is at index 1.
    xmax=50; xmin=-50; ymax=50; ymin=-50; %Environment size.
    dxv=[0.35; 0.35];  %Give the animal starting momentum so it won't wander in loops at start.
    for t=1:fullT
        dxv(:,t+1)=((1-mom)*speed*randn(2,1))/dt + mom*dxv(:,t); %Movement in x and y direction (2D vector) need to divide by dt for true velocity.
        xv(:,t+1) = xv(:,t) + dt*dxv(:,t+1);
        %BARRIERS
        if xv(1,t+1)<xmax && xv(1,t+1)>xmin;
        else; dxv(1,t+1)=-revs*dxv(1,t+1); xv(1,t+1)=xv(1,t)+dt*dxv(1,t+1); end; %If hits barrier, give reverse momentum term.
        if xv(2,t+1)<ymax && xv(2,t+1)>ymin;
        else; dxv(2,t+1)=-revs*dxv(2,t+1); xv(2,t+1)=xv(2,t)+dt*dxv(2,t+1); end 
    end
    pos_x=xv(1,:);
    pos_y=xv(2,:);
    espeed=zeros(1,fullT+1);
    for t=1:fullT
        espeed(1,t)=sqrt(dxv(1,t+1)^2+dxv(2,t+1)^2);
        sumspeed=sumspeed+espeed(1,t);
    end
    sumspeed*dt %Print total distance covered.
    fullT*dt     %Print total time run.
    avespeed=sumspeed/fullT %Print average speed.
    

elseif datainput==4;  %Circle without stopping points.
    dxv=[0.35; 0.35];
    dt=0.02; sumspeed=0;
    xmax=60; xmin=-60; ymax=60; ymin=-60; %Environment size for circular track.
    radius=47.5;
    for t=1:fullT
        if t==1 xv=[radius; 0]; phi(:,1)=(HHst*xv(:,1)>0).*(HHst*xv(:,1)); 
        end; 
        cirangle=atan2(xv(2,t),xv(1,t));
        dxv(1,t+1)=(1.0+0.0*rand)*sin(cirangle-1/96)/dt; %/dt; %1/4 larger inward vs. 1/ 64 larger outward.  Use 1/32 to get larger inward movement?
        dxv(2,t+1)=-(1.0+0.0*rand)*cos(cirangle-1/96)/dt; %/dt; %For 20*sin -1/5 is best.  for 10*sin
        xv(:,t+1) = xv(:,t) + dt*dxv(:,t+1) ;
        speed=sqrt(dxv(1,t+1)^2+dxv(2,t+1)^2);
        sumspeed=sumspeed+speed;
    end
    sumspeed*dt  %Print total distance covered.
    fullT*dt  %Print total time run.
avespeed=sumspeed/fullT  %Print average speed.
pos_x=xv(1,:);
pos_y=xv(2,:);
elseif datainput==5; %Circle with stopping points.
    dxv=[0; -0.35];
    reward=zeros(1,fullT);
    espeed=zeros(1,fullT+1);
    dt=0.02; sumspeed=0;
    xmax=60; xmin=-60; ymax=60; ymin=-60; %Environment size for circular track.
    radius=47.5;
    stopangle=0; 
    ccstop=1;
    for t=1:fullT
        if t==1 xv=[radius; 0]; phi(:,1)=(HHst*xv(:,1)>0).*(HHst*xv(:,1)); 
        end; 
        cirangle=atan2(xv(2,t),xv(1,t));
        if t>1
        cirangleOLD=atan2(xv(2,t-1),xv(1,t-1));
        dangle=cirangleOLD-cirangle;
        if abs(dangle)<pi/2
        stopangle=stopangle+dangle;
        end
        end
        if xv(2,t)<12 && xv(2,t)>2 && xv(1,t)>20
            if dur==0; dur=round(rand*1000); dur
            end;
            ccstop=ccstop+1;
            dxv(1,t+1)=(0.04+0.02*rand)*sin(cirangle-1/96)/dt; %This causes Slower movement at reward locations
            dxv(2,t+1)=-(0.04+0.02*rand)*cos(cirangle-1/96)/dt;
            [ccstop t+1 xv(1,t) xv(2,t)];
            preRewStop=50;
            if ccstop>preRewStop && ccstop<(preRewStop+dur); %Cause stopping for a period of time.
                dxv(1,t+1)= 0 + 0.02*rand*0.1/dt;   % +  0.02*rand*sin(cirangle-1/96)/dt;
                dxv(2,t+1)= 0 + 0.02*rand*0.1/dt;   %*cos(cirangle-1/96)/dt;
                reward(1,t)=1;
                stopangle=0;
                %ccstop=1;
            end
            xv(:,t+1) = xv(:,t) + dt*dxv(:,t+1) ; %This updates with slow movement
        else
        dxv(1,t+1)=(1.0+0.02*rand)*sin(cirangle-1/96)/dt; %/dt; %1/4 larger inward vs. 1/ 64 larger outward.  Use 1/32 to get larger inward movement?
        dxv(2,t+1)=-(1.0+0.02*rand)*cos(cirangle-1/96)/dt; %/dt; %For 20*sin -1/5 is best.  for 10*sin
        xv(:,t+1) = xv(:,t) + dt*dxv(:,t+1) ;
        dur=0; ccstop=1; %Reset counters for duration of stopping period.
        end
        espeed(1,t+1)=sqrt(dxv(1,t+1)^2+dxv(2,t+1)^2);
        sumspeed=sumspeed+espeed(1,t+1);
    end
    sumspeed*dt  %Print total distance covered.
    fullT*dt  %Print total time run.
avespeed=sumspeed/fullT  %Print average speed.
pos_x=xv(1,:);
pos_y=xv(2,:);
end %End all datainput

if datainput~=7 && datainput~=8;
    Stim=[0 0 0 0 0;
          0 0 0 0 0;
          0 0 0 0 0;
          0 0 0 0 0;
          0 0 0 0 0];
    StimVect=[0 0 0 0 0;
              0 0 0 0 0;
              0 0 0 0 0;
              0 0 0 0 0;
              0 0 0 0 0];
end

for et=3:fullT;
    espeed(1,et)=sqrt((xv(1,et)-xv(1,et-1))^2 + (xv(2,et)-xv(2,et-1))^2)/dt;
end

%SOME GRAPHICS PARAMETERS
plotH=0; %Sets whether head direction vectors are plotted at end.
gy=0.7; %Sets trajectory color to gray.
figsize=200;
plotGrids=0;
plotCells=0;
plotEnc=1;
plotRet=1;

if plotEnc==1
axes(ax111);
tt=[1:fullT];
plot(tt,xv(2,1:fullT),tt,espeed(1,1:fullT));
xlim([1 fullT]);
ylim([ymin-10 ymax+10]);

%PLOT ACTUAL TRAJECTORY IN GRAY
axes(ax14b);
gy=0.7;
set(gcf,'DefaultAxesColorOrder',[gy gy gy])
plot(pos_x(minT:minT+fullT), pos_y(minT:minT+fullT),'LineWidth',2); hold on;
plot(Stim(1,1), Stim(1,2),'+','MarkerSize',10,'LineWidth',2); hold on;
plot(Stim(2,1), Stim(2,2),'*','MarkerSize',10,'LineWidth',2); hold on;
plot(Stim(3,1), Stim(3,2),'o','MarkerSize',10,'LineWidth',2); hold on;
plot(Stim(4,1), Stim(4,2),'s','MarkerSize',10,'LineWidth',2); hold on;
plot(Stim(5,1), Stim(5,2),'^','MarkerSize',10,'LineWidth',2); hold on;
xlim([xmin-10 xmax+10]);
ylim([ymin-10 ymax+10]);
end

%******INITIAL CELL CREATION************************
%CREATE GRID CELLS FOR INITIAL PLACE CELL CREATION

T=zeros(1,fullT);  T(1,1)=1;
iphi=[0; 0; 0]; %Vector of initial phases of dendrites (in radians) %iphi=[pi -0.5*pi -0.5*pi];
B=0.00385; %Scaling factor for effect of head direction on membrane oscillation frequency.
%strengHD=0; %Strength of head direction input during retrieval.
numF=4; %Number of different frequencies (usually 3 in most simulations)
baseF=0; %Baseline frequency for range. (usually 0 in most simulations using cos(phi-phis)
multFF=1; %Addition for going through frequencies.
numPhi=5; %Number of different starting phases in each dimension (default=5 for most simulations)
%numP=20; %Number of rows (and columns) of place cells.
thresh=0.3; %Threshold for grid cells.
pOff=0.0; %Grid cell phase offset.
stdlim=2.7; %2.7 best for circle %Determines size of accepted place fields.
xspike=zeros(2,numF,numPhi^2,fullT); %Location of grid cellspiking
tspike=zeros(numF*numPhi^2,fullT); %Spiking of grid cells during creation.
tspike2=zeros(1,numF*numPhi^2,fullT); %Spiking of grid cells during encoding.
timep(1,numF*numPhi^2).spike=[];
allPlotPlaceX=zeros(1,1); %fullT
allPlotPlaceY=zeros(1,1);
sumPlace=zeros(1,fullT);

if plotGrids==1
figure('name','grid creation','Numbertitle','off','Position',[250 50 numPhi*150 numPhi*150]); ax15=gca;
end

%CREATION OF GRID CELLS AND PLACE CELLS.

useDXV=1;
%Cycle through frequencies and phases for creating grid cells.
for ff=1:numF; %GRID CELL INITIAL CREATION.
  f=baseF+multFF*ff;  
  for pp=1:numPhi;
    for p2=1:numPhi;
        
phi=zeros(HDnum,fullT); %Vector for phase of dendrites (in radians).
phis=zeros(HDnum/2,fullT);  %Same soma phase in all parts of vector, but easier to compute with this.

    ix=[xv(1,1)+((pp-1)+(ff-1)+pOff)*(sqrt(3)/(f*B*(numPhi))); xv(2,1)+((p2-1)+(ff-1)+pOff)*(sqrt(3)/(f*B*(numPhi)))]; %f*G=H=sqrt(3)/2*B. G=sqrt(3)/2*f*B
    hx=((HHst*ix)>0).*(HHst*(ix));
    %LH=length(h)/2; MH=length(h);
       if useDXV==1; %phis(:,1)=2*pi*f*dt; 
           phi(:,1)=2*pi*f*(B*hx); end;

         for t=2:fullT;
           if useDXV==1; %phis(:,t)=phis(:,t-1)+2*pi*f*dt; 
               phi(:,t)=phi(:,t-1)+2*pi*f*dt*(1+B*((HHst*dxv(:,t))>=0).*(HHst*dxv(:,t))); 
           %elseif useDXV==0; phis(:,t)=2*pi*f*t; phi(:,t)=2*pi*f*(t+B*HHst*(xv(:,t)-ix));
           end; %Creation of phase representation.
         end

LH=length(h)/2; MH=length(h);
g=prod(cos((phi(1:LH,:)-phi(LH+1:MH,:))));

spiket=g>thresh; %Spiking of grid cell.

tspike((ff-1)*numPhi^2+(pp-1)*numPhi+p2,1:fullT)=spiket; %Store spiking of all grid cells.

if plotGrids==1 %OPTIONAL GRID CELL FIGURE
axes(ax15); 
axes(ax15);
plotSpikeX=xv(1,1:fullT).*spiket;
plotSpikeY=xv(2,1:fullT).*spiket;
subplot(numPhi,numPhi,(p2-1)*numPhi+pp,'replace'); ax15=gca;
plot(plotSpikeX,plotSpikeY,'r.');
xlim([xmin-1 xmax+1])
ylim([ymin-1 ymax+1])
set(gca,'XTickLabel',{' '})
set(gca,'YTickLabel',{' '})
drawnow
end %End if plotGrids==1

toc;
        end %End for p2=1:numPhi
    end %End for pp=1:numPhi
end %End for f=1:numF

%INITIAL PLACE CELLS CREATION - place cells.
%numP=6; Number set with numP above.
placeSpike=zeros(numP^2,fullT);
if plotCells==1
figure('name','place cells creation','Numbertitle','off','Position',[50 50 800 800]); ax16=gca;
set(gcf,'DefaultAxesColorOrder',[gy gy gy])
set(gca,'XTickLabel',{' '})
set(gca,'YTickLabel',{' '})
drawnow
end %End if plotCells==1

WPG=zeros(numP^2,numF*(numPhi^2)); %WPG matrix will store associations between grid cells and place cells.
numGtoP=3;
nP=1;
nnnn=1;
while nP<(numP^2+1)
    nnnn=nnnn+1;
    if nnnn>50000; break; end;
    gp=1;
    pPhi=ones(gp,1);
for ppp=1:numF*numPhi^2;  %SELECT THREE GRID CELLS AT RANDOM WITHOUT OVERLAP.
    pRand=rand>0.9; %Needs to be 0.9 or larger.  0.8 made it VERY slow - too many.
    if pRand==1; pPhi(gp,1)=ppp; gp=gp+1; end; %Sequentially activates elements.
    if gp==numGtoP+1; break; end; %ppp=1:numF*numPhi^2; end;
end

pspike=prod(tspike(pPhi,1:fullT));  %Multiplies together random set of grid cells.
plotPlaceX=xv(1,1:fullT).*pspike;
plotPlaceY=xv(2,1:fullT).*pspike;

if length(find(pspike))>3 && gp==numGtoP+1; %Only take place cells with at least three spikes from three grid inputs gp.

    stX=std(plotPlaceX(1,find(pspike)));
    stY=std(plotPlaceY(1,find(pspike)));
    
    if stX<stdlim && stY<stdlim && (stX>0.1 || stY>0.1); %Selects only place cells with three grid inputs and >3 spikes
        
placeSpike(nP,:)=pspike;

if plotCells==1
%Plot place spike.
axes(ax16);
subplot(numP,numP,nP,'replace'); ax16=gca;  %Plots the place cell spike location plots

plot(plotPlaceX,plotPlaceY,'r.');
xlim([xmin-1 xmax+1])
ylim([ymin-1 ymax+1])
set(gca,'XTickLabel',{' '})
set(gca,'YTickLabel',{' '})
end %End plotCells==1

%Keep track of times of spikes of selected place cells.
sumPlace=sumPlace+pspike;

WPG(nP,pPhi)=1; %For a successful place cell, adds connectivity to WPG for input from grids.
    nP=nP+1;
    
    end %End std loop
end %End find(pspike) loop
end %End nP loop.

nnnn
toc;

%Plot the sum of place cell firing.
plotAllPlace=0;

if plotAllPlace==1;
    %Always plot the coverage of all place cells
figure('name','all place cells','Numbertitle','off','Position',[50 500 200 200]); ax17=gca;
set(gcf,'DefaultAxesColorOrder',[gy gy gy])
axes(ax17);
plot(xv(1,1:fullT).*(sumPlace>0.5),xv(2,1:fullT).*(sumPlace>0.5),'ko'); %Plots the sum of all place cell firing to see coverage.
xlim([xmin-10 xmax+10]);
ylim([ymin-10 ymax+10]);
drawnow
end

%********ENCODING OF TRAJECTORY WITH PLACE, HD AND GRID************************
eeeT=fullT; %100; %fullT;
tspike2=zeros(numF*numPhi^2,eeeT);
placeSpike2=zeros(numP^2,eeeT);

arctspike2=zeros(arcnum,eeeT);
placeArcSpike2=zeros(numP^2+arcnum,eeeT);
%WHP=zeros(HDnum,numP^2+1); 
WHP=zeros(HDnum,numP^2+arcnum); %Creates matrix from place cells and arc cells to HD cells.
WSP=zeros(length(Stim(:,1)),numP^2+arcnum);
WHA=zeros(1,arcnum);
WPA=zeros(HDnum,arcnum);
WPAs=zeros(HDnum,arcnum);

pephi=zeros(length(h),eeeT);
plotephi=zeros(length(h)/2,eeeT);
plothd=zeros(length(h)/2,eeeT);

distStim=zeros(length(Stim(:,1)),eeeT); 

%Grid cell phase arrays

arcphis=zeros(arcnum,eeeT);
arcphi=zeros(arcnum,eeeT);

useDXV=1;
%Cycle through frequencies and phases for creating grid cells.
for ff=1:numF; %GRID CELLS DURING ENCODING (USING SAME AS DURING CREATION.
    f=baseF+multFF*ff;  
    for pp=1:numPhi;
        for p2=1:numPhi;
            
ephis=zeros(length(h),eeeT);
ephi=zeros(length(h),eeeT);
ephiOLD=zeros(length(h),eeeT);
ephisOLD=zeros(length(h),eeeT);

        ix=[xv(1,1)+((pp-1)+(ff-1)+pOff)*(sqrt(3)/(f*B*(numPhi))); xv(2,1)+((p2-1)+(ff-1)+pOff)*(sqrt(3)/(f*B*(numPhi)))]; %f*G=H=sqrt(3)/2*B. G=sqrt(3)/2*f*B
        hx=((HHst*ix)>0).*(HHst*(ix));
        %LH=length(h)/2; MH=length(h);
        if useDXV==1; ephi(:,1)=2*pi*f*(B*hx); end;
        if useDXV==1 && ff==1 && pp==1 && p2==1; pephi(:,1)=2*pi*f*(B*hx); end;
        
for t=2:eeeT;
if useDXV==1; ephis(:,t)=ephis(:,t-1)+2*pi*f*dt; 
    ephi(:,t)=ephi(:,t-1)+2*pi*f*dt*(1+B*((HHst*dxv(:,t))>=0).*(HHst*dxv(:,t)));
    if ff==1 && pp==1 && p2==1;     
      pephi(:,t)=pephi(:,t-1)+2*pi*f*dt*(1+B*((HHst*dxv(:,t))>=0).*(HHst*dxv(:,t))); 
    end
elseif useDXV==0; ephis(:,t)=2*pi*f*t; ephi(:,t)=2*pi*f*(t+B*HHst*(xv(:,t)-ix));
end; %Creation of phase representation.
end

LH=length(h)/2; MH=length(h);
g=prod(cos((ephi(1:LH,:)-ephi(LH+1:MH,:))));

spiket=g>thresh; %Spiking of grid cell.

tspike2((ff-1)*numPhi^2+(pp-1)*numPhi+p2,1:fullT)=spiket; %Store spiking of all grid cells.

        end %End for p2=1:numPhi
    end %End for pp=1:numPhi
end %End for f=1:numF

for tt=2:fullT;
 plotephi(:,tt)=pephi(1:LH,tt)-pephi(LH+1:MH,tt);
end
     

farc=0.6; %0.6; %0.6*500/fullT;
arcthresh=0.9999;


for t=1:eeeT;

%PLACE CELLS DURING ENCODING
placeSpike2(:,t)=((WPG*tspike2(:,t))./numGtoP)>0.9; %Spiking of place cells during encoding.

speed(1,t)=sqrt(dxv(1,t)^2+dxv(2,t)^2);


%CREATE ARC LENGTH CELLS

if t==1;  
    for arcn=1:arcnum;
        arcphi(arcn,1)=2*pi*arcn/arcnum;
    end
elseif t>1;
    arcspeed=sqrt(dxv(1,t)^2+dxv(2,t)^2);
    for arcn=1:arcnum;
         if useDXV==1; 
             arcphis(arcn,t)=arcphis(arcn,t-1)+2*pi*farc*dt; 
             arcphi(arcn,t)=arcphi(arcn,t-1)+2*pi*farc*dt*(1+B*arcspeed); 
         end; %Creation of phase representation.
         arctspike2(arcn,t)=cos(arcphis(arcn,t)-arcphi(arcn,t))>arcthresh;
     end %End arcn
     
end %End if t==1, t>1
 
%HD CELLS DURING ENCODING
hd(:,t)=((HHst*dxv(:,t))>=0).*(HHst*dxv(:,t)); %Add the head direction cell activity during encoding.

hd(:,t);

    plothd(:,t)=hd(1:LH,t)-hd(LH+1:MH,t);

%PLACE TO HD CELL
placeArcSpike2(1:numP^2,t)=placeSpike2(:,t);
placeArcSpike2(numP^2+1:numP^2+arcnum,t)=arctspike2(:,t);
%placeArcSpike2=[placeSpike2; arctspike2];
norm=sum(placeArcSpike2(:,t));
if norm~=0;
    for xx=1:numP^2+arcnum; %Be sure to include arc cells here!!
        if placeArcSpike2(xx,t)==1;
WHP(:,xx)=(WHP(:,xx)+hd(:,t).*placeArcSpike2(xx,t))/2; %Create connections from place to head direction.
%WHP(:,xx)=(WHP(:,xx)+hd(:,t).*placeSpike2(xx,t)./norm)/2; %Create connections from place to head direction.
        end  %This procedure takes elements with presynaptic input and averages the existing and new weight.
    end
end


%PLACE CELLS TO STIM
StimLim=12;
for ss=1:length(Stim(:,1));
  distStim(ss,t)=(xv(1,t)-Stim(ss,1))^2 + (xv(2,t)-Stim(ss,2))^2;
  if distStim(ss,t) < (StimLim^2)
    for xx=1:numP^2 + arcnum;
        if placeArcSpike2(xx,t)==1;
    newWSP=(StimVect(:,ss).*placeArcSpike2(xx,t)); %placeArcSpike2
    WSP(:,xx)= (WSP(:,xx) + newWSP)/2;
        end
    end
  end
end


%ARC CELL TO HD CELL
normA=sum(arctspike2(:,t));
if normA~=0;
    for xx=1:arcnum;
        if arctspike2(xx,t)==1 %&& max(arctspike2(xx,t).*arctspike2(xx,t-1))~=1;
            WHA(:,xx)=(WHA(:,xx)+speed(:,t).*arctspike2(xx,t))/2; 
        end%Dividing by /2 gives an average of current
        if arctspike2(xx,t)==1 && max(arctspike2(xx,t).*arctspike2(xx,t-1))~=1;
            WPA(:,xx)=pephi(:,t).*arctspike2(xx,t); %removed WPA(:,xx)+  %%divide by HDnum?
            WPAs(:,xx)=ephis(:,t).*arctspike2(xx,t);
        end
    end
end

end %End for t=1:eeeT


if plotGrids==1
figure('name','grid during encoding','Numbertitle','off','Position',[250 50 numPhi*150 numPhi*150]); ax15b=gca;
axes(ax15b); %Grid cell figure

%Plot Grid cells during encoding.
for ff=1:numF;
    f=baseF+multFF*ff;
    for pp=1:numPhi;
        for p2=1:numPhi;
plotSpikeX=xv(1,1:eeeT).*tspike2((ff-1)*numPhi^2+(pp-1)*numPhi+p2,1:eeeT);
plotSpikeY=xv(2,1:eeeT).*tspike2((ff-1)*numPhi^2+(pp-1)*numPhi+p2,1:eeeT);
axes(ax15b);
%subplot(numF*numPhi,numPhi,(f-1)*numPhi^2+(pp-1)*numPhi+p2); ax15b=gca;
subplot(numPhi,numPhi,(pp-1)*numPhi+p2,'replace'); ax15b=gca;
plot(plotSpikeX,plotSpikeY,'r.');
xlim([xmin-1 xmax+1])
ylim([ymin-1 ymax+1])
set(gca,'XTickLabel',{' '})
set(gca,'YTickLabel',{' '})
drawnow
       end %End for p2=1:numPhi
    end %End for pp=1:numPhi
end %End for f=1:numF
end %End if plotGrids==1

%Plot Arc cells during encoding.
plotArcs=0;
if plotArcs==1
figure('name','arc cells encoding','Numbertitle','off','Position',[50 50 1100 60]); ax15bb=gca;
axes(ax15bb); %Arc cell figure

for arcn=1:arcnum;
    plotArcSpikeX=xv(1,1:eeeT).*arctspike2(arcn,1:eeeT);
    plotArcSpikeY=xv(2,1:eeeT).*arctspike2(arcn,1:eeeT);
axes(ax15bb);
subplot(1,arcnum,arcn); ax15bb=gca;
plot(plotArcSpikeX,plotArcSpikeY,'r.');
xlim([xmin-1 xmax+1])
ylim([ymin-1 ymax+1])
set(gca,'XTickLabel',{' '})
set(gca,'YTickLabel',{' '})
drawnow%
end  %end arcn.
end %if plotArc==1

if plotCells==1
figure('name','place cells Encoding','Numbertitle','off','Position',[800 50 numP*100 numP*100]); ax16bbb=gca;

%Plot Place Cells during encoding.
for pp=1:numP^2;
plotPlaceX=xv(1,1:eeeT).*placeSpike2(pp,:);
plotPlaceY=xv(2,1:eeeT).*placeSpike2(pp,:);
axes(ax16bbb);
subplot(numP,numP,pp,'replace'); ax16bbb=gca;
plot(plotPlaceX,plotPlaceY,'r.');
xlim([xmin-1 xmax+1])
ylim([ymin-1 ymax+1])
set(gca,'XTickLabel',{' '})
set(gca,'YTickLabel',{' '})
end %End for pp
end %End if plotCells==1

if plotEnc==1;
    
%PLOT PLACE CELLS VS TIME DURING ENCODING spike times with colors.
axes(ax16b);
subplot(numPsizesTot*numPtestsTot,1,(numPsizes-1)*numPtestsTot+numPtests,'replace'); ax16b=gca;

numSING=numPlaceToPlot;
for sing=1:numSING;
    fpsP=find(placeSpike2(sing,:));
    for mm=1:length(fpsP);
        if colorSpikes==1;
            plot([fpsP(mm) fpsP(mm)],[sing sing-0.5],'Color',[rc(sing,:)],'LineWidth',1); hold on;
        elseif colorSpikes==0;
            plot([fpsP(mm) fpsP(mm)],[sing sing-0.5],'Color','k','LineWidth',1); hold on;
        end
    end
end
if colorSpikes==1;
    set(gca,'Color','k');
elseif colorSpikes==0
    set(gca,'Color','w');
end
xlim([1 eeeT]);
ylim([0 numSING+0.5]);
drawnow

%if plotEnc==1;

%PLOT HD cells during encoding
axes(ax18);

for m=1:20:eeeT-1;   %Plot the head direction cell vectors. %If want to plot those during encoding start m=1:
     denom=cos(h(1,1))*sin(h(2,1))-sin(h(1,1))*cos(h(2,1));
     hdx=((plothd(1,m+1))*sin(h(2,1))-(plothd(2,m+1))*sin(h(1,1))); %./(denom*wd*B);
     hdy=((plothd(1,m+1))*-cos(h(2,1))+(plothd(2,m+1))*cos(h(1,1))); %./(denom*wd*B);
    
      AA=[xv(1,m) xv(1,m)+hdx];  
      BB=[xv(2,m) xv(2,m)+hdy];
      
      CC=[xv(1,m) xv(1,m)+dxv(1,m+1)*dt];
      DD=[xv(2,m) xv(2,m)+dxv(2,m+1)*dt];
      if sum(AA)~=0 && sum(BB)~=0;
      plot(AA,BB,'k-',CC,DD,'r-','LineWidth',1); %AA,BB,'k-',
      hold on
      end
end
xlim([xmin-10 xmax+10]);
ylim([ymin-10 ymax+10]);

%PLOT HD ANGLE VERSUS TIME
axes(ax113);
tt=[1:fullT];
plot(tt,plothd(1,:));
xlim([1 fullT]);
ylim([-20 20]);

%PLOT HD CELLS VERSUS TIME DURING ENCODING
axes(ax18b)

fpsHDeee.cell=[];
spikehd=[];
spikehd=zeros(length(h),length(hd));
for ht=1:20:length(hd); 
    maxhd=mean(max(hd));
    spikehd(:,ht)=rand(length(h),1)>(maxhd*ones(length(h),1)-hd(:,ht))/maxhd; 
end
numSING=length(h);
for sing=1:numSING;
    fpsE=find(spikehd(sing,:));
    fpsHDeee(sing).cell=find(spikehd(sing,:));
    for mm=1:length(fpsE);
%plot([fps(mm) fps(mm)],[sing sing-0.5],'Color',[rc(sing,:)],'LineWidth',1); hold on;
        if colorSpikes==1;
            plot([fpsE(mm) fpsE(mm)],[sing sing-0.5],'Color',[rc(sing,:)],'LineWidth',1.8); hold on;
        elseif colorSpikes==0;
            plot([fpsE(mm) fpsE(mm)],[sing sing-0.5],'Color','k','LineWidth',1); hold on;
        end
    end
end
if colorSpikes==1;
    set(gca,'Color','k');
elseif colorSpikes==0
    set(gca,'Color','w');
end
xlim([1 eeeT]);
ylim([0 numSING+0.5]);
drawnow

end %End plotEnc==1

toc;

%********RETRIEVAL OF TRAJECTORY WITH PLACE, HD AND GRID************************
%Add RETRIEVAL period.

rrrT=round(fullT); %0.2*fullT; %eeeT;
hd3=zeros(length(h),rrrT);
stim3=zeros(length(Stim(:,1)),rrrT);
hd3(:,1:5)=hd(:,5:9); %Set first step of retrieval hd to same as encoding hd.
tspike3=zeros(numF*numPhi^2,rrrT); 
tspike3(:,1:5)=tspike2(:,5:9); %Set initial grid cell spiking to be same as during encoding
%Grid cell spiking initial conditions don't matter. Can set to tspike2(:,21:22) and still function.  But does initial phase set by ix matter?
placeSpike3=zeros(numP^2,rrrT);
placeSpike3(:,1:5)=placeSpike2(:,5:9); %set initial place cell to be same as during encoding.
arctspike3=zeros(arcnum,rrrT);
placeArcSpike3=zeros(numP^2+arcnum,rrrT);
arcphis3=zeros(arcnum,rrrT);
arcphi3=zeros(arcnum,rrrT);
arcspeed3=zeros(1,rrrT);

rphis=zeros(length(h),numF*numPhi^2); 
rphi=zeros(length(h),numF*numPhi^2); 
aRphi=zeros(HDnum,rrrT);
aRphis=zeros(round(HDnum),rrrT);
plotrphi=zeros(length(h)/2,rrrT);
plotiG=zeros(length(h),1);

rphiOLD=zeros(length(h),numF*numPhi^2);
rphisOLD=zeros(length(h)/2,numF*numPhi^2);

if useDXV==1; %Set initial conditions for grid cells.
for ff=1:numF; 
    f=baseF+multFF*ff;
    for pp=1:numPhi;
        for p2=1:numPhi;
            cG=(ff-1)*numPhi^2+(pp-1)*numPhi+p2; %Number of current grid cell.
        ix=[xv(1,1)+((pp-1)+(ff-1)+pOff)*ceil(sqrt(3)/(f*B*(numPhi))); xv(2,1)+((p2-1)+(ff-1)+pOff)*ceil(sqrt(3)/(f*B*(numPhi)))]; %f*G=H=sqrt(3)/2*B. G=sqrt(3)/2*f*B
           iG=((HHst*ix)>0).*(HHst*ix);
           for nG=1:length(h);
             %rphis(nG,cG)=2*pi*f*dt; 
             rphi(nG,cG)=2*pi*f*(B*iG(nG)); 
           end
           if cG==1
               plotiG=(HHst*(ix-xv(:,1))>0).*(HHst*(ix-xv(:,1)));
           end
           end %End for p2=1:numPhi
    end %End for pp=1:numPhi
end %End for f=1:numF
end %End if useDXV.

%SET INITIAL ARC CELL PHASE
for arcn=1:arcnum;
    arcphi3(arcn,1)=2*pi*arcn/arcnum;
    arctspike3(arcn,1)=cos(0-arcphi3(arcn,1))>arcthresh;
end

%RETRIEVAL ACTIVITY OF GRID CELLS, PLACE CELLS AND HD CELLS
for t=2:rrrT;

    
    %GRID CELLS
    for ff=1:numF; 
    f=baseF+multFF*ff;
    for pp=1:numPhi;
        for p2=1:numPhi;
            cG=(ff-1)*numPhi^2+(pp-1)*numPhi+p2; %Number of current grid cell.
        ix=[xv(1,1)+((pp-1)+(ff-1)+pOff)*ceil(sqrt(3)/(f*B*(numPhi))); xv(2,1)+((p2-1)+(ff-1)+pOff)*ceil(sqrt(3)/(f*B*(numPhi)))]; %f*G=H=sqrt(3)/2*B. G=sqrt(3)/2*f*B

    %RETRIEVE GRID PHASE FROM ARC CELLS
    arcToRphi=0; %Keep at zero as it does not seem to work.
    if arcToRphi==1;
    rAphi(:,t)=WPA*arctspike3(:,t-1); %/HDnum; %sqrt(retdxv(1,1)^2+retdxv(2,1)^2);
    %rAphis(:,t)=WPAs*arctspike3(:,t-1); %/HDnum; 
        hx=((HHst*ix)>0).*(HHst*(ix));
    end
    
%IF ARCTOPLACE==1, THEN UPDATE GRID PHASE WITH MATRIX WPA.
%arcToRphi=1;
if t>2 && arcToRphi==1 && max(arctspike3(:,t-1))==1 && max(arctspike3(:,t-1).*arctspike3(:,t-2))~=1; 
    
        if useDXV==1; %rphis(:,cG)=rAphis(:,t); 
            rphi(:,cG)=rAphi(:,t); end; %+2*pi*f*(B*hx); end;
end
        
        %UPDATE grid phase based on head direction
if useDXV==1; 
    rphiOLD=rphi;
    rphi(:,cG)=rphiOLD(:,cG)+2*pi*f*dt*(1+B*strengHD*hd3(:,t));  %Use head direction vector hd3 retrieved from WHP.

end



LH=length(h)/2; MH=length(h);
tspike3(cG,t)=prod(cos(rphi(1:LH,cG)-rphi(LH+1:MH,cG)))>thresh; 

if cG==1;
    plotrphi(:,t)=rphi(1:LH,cG)-rphi(LH+1:MH,cG);
    plothd3(:,t)=hd3(1:LH,t)-hd3(LH+1:MH,t);
end
       end %End for p2=1:numPhi
    end %End for pp=1:numPhi
end %End for f=1:numF

%PLACE CELLS DURING RETRIEVAL
placeSpike3(:,t)=((WPG*tspike3(:,t))./numGtoP)>0.9; %Create place cell activity during encoding.
    

if t>1;
    arcspeed3(1,t)=WHA*arctspike3(:,t-1); %sqrt(retdxv(1,1)^2+retdxv(2,1)^2);
    for arcn=1:arcnum;
         if useDXV==1; 
             arcphis3(arcn,t)=arcphis3(arcn,t-1)+2*pi*farc*dt; 
             arcphi3(arcn,t)=arcphi3(arcn,t-1)+2*pi*farc*dt*(1+B*arcspeed);
         end; %Creation of phase representation.
         arctspike3(arcn,t)=cos(arcphis3(arcn,t)-arcphi3(arcn,t))>arcthresh;
     end %End arcn
end %End if t==1, t>1

%COMBINE PLACE CELLS AND ARC LENGTH CELLS
%placeArcSpike3(:,t)=[placeSpike3(:,t); arctspike3(3,t)];
placeArcSpike3(:,t)=[placeSpike3(:,t); arctspike3(:,t)];

%HD CELLS
if max(placeArcSpike3(:,t))==0
    hd3(:,t+1)=hd3(:,t);
else 
    norm=sum(placeArcSpike3(:,t));
hd3(:,t+1)=WHP*placeArcSpike3(:,t)/norm;
%RETRIEVAL OF STIM
stim3(:,t+1)=WSP*placeArcSpike3(:,t)/norm;
end

end %End for t=1:rrrT

%PLOT INVERSE TRANSFORM TO SEE WHERE RAT "THINKS" IT IS LOCATED.
axes(ax14); %Go back to x y plot and plot inverse.
%subplot(numPsizesTot, numPtestsTot, (numPsizes-1)*numPtestsTot+numPtests); ax14=gca; %This does rows of tests.
subplot(numPsizesTot,numPtestsTot,(numPsizes-1)*numPtestsTot+numPtests,'replace'); ax14=gca; %This does one column of all.

%PLOT ACTUAL TRAJECTORY IN GRAY
set(gcf,'DefaultAxesColorOrder',[gy gy gy])

plot(pos_x(minT:minT+fullT), pos_y(minT:minT+fullT),'LineWidth',2);
hold on


%PLOT RETRIEVED TRAJECTORY IN BLACK
set(gcf,'DefaultAxesColorOrder',[0 0 0])

 wd=(baseF+multFF*1)*2*pi;
 
 %Calculate inverse transform of grid cell phase (plotrphi(1,:) and plotrphi(2,:)
    denom=cos(h(1,1))*sin(h(2,1))-sin(h(1,1))*cos(h(2,1));
    xphE=((plotephi(1,:))*sin(h(2,1))-(plotephi(2,:))*sin(h(1,1)))./(denom*wd*B);
    yphE=((plotephi(1,:))*-cos(h(2,1))+(plotephi(2,:))*cos(h(1,1)))./(denom*wd*B);
    xph=((plotrphi(1,:))*sin(h(2,1))-(plotrphi(2,:))*sin(h(1,1)))./(denom*wd*B);
    yph=((plotrphi(1,:))*-cos(h(2,1))+(plotrphi(2,:))*cos(h(1,1)))./(denom*wd*B);
    plot(xphE(1,2:rrrT),yphE(1,2:rrrT));
    plot(xph(1,2:rrrT),yph(1,2:rrrT),'k','LineWidth',2); %xlim([xmin+30 xmax+30]); ylim([ymin-60 ymax-60]); 


    
%Plot STIMULUS RETRIEVAL

plot(Stim(1,1), Stim(1,2),'+','MarkerSize',10,'LineWidth',2); hold on; %was blue 'bo'
plot(Stim(2,1), Stim(2,2),'*','MarkerSize',10,'LineWidth',2); hold on;
plot(Stim(3,1), Stim(3,2),'o','MarkerSize',10,'LineWidth',2); hold on;
plot(Stim(4,1), Stim(4,2),'s','MarkerSize',10,'LineWidth',2); hold on; %was blue 'bo'
plot(Stim(5,1), Stim(5,2),'^','MarkerSize',10,'LineWidth',2);
%plot(Stim(3,1), Stim(3,2),'o','MarkerSize',10,'LineWidth',2);

for tt=5:rrrT
    sth=0.95;
    if stim3(1,tt)>sth && max(stim3(1,tt-2:tt-1))<sth;
        plot(xph(1,tt)-2,yph(1,tt),'k+','MarkerSize',12,'LineWidth',2); hold on;
    elseif stim3(2,tt)>sth && max(stim3(2,tt-2:tt-1))<sth;
        plot(xph(1,tt)-2,yph(1,tt),'k*','MarkerSize',12,'LineWidth',2); hold on;
    elseif stim3(3,tt)>sth && max(stim3(3,tt-2:tt-1))<sth;
        plot(xph(1,tt)-2,yph(1,tt),'ko','MarkerSize',12,'LineWidth',2); hold on; %was red 'r+'
    elseif stim3(4,tt)>sth && max(stim3(4,tt-2:tt-1))<sth;
        plot(xph(1,tt)-2,yph(1,tt),'ks','MarkerSize',12,'LineWidth',2); hold on;
    elseif stim3(5,tt)>sth && max(stim3(5,tt-2:tt-1))<sth;
        plot(xph(1,tt)-2,yph(1,tt),'k^','MarkerSize',12,'LineWidth',2); 
    %elseif stim3(3,tt)>sth && max(stim3(3,tt-2:tt-1))<sth;
    %    plot(xph(1,tt)-2,yph(1,tt),'ko','MarkerSize',12,'LineWidth',2);
    end
end

    
xlim([xmin-10 xmax+10])
ylim([ymin-10 ymax+10])
if datainput==1
xlim([-85 -15]);
ylim([-25 65]);
end
%xlim([-500 500])
%ylim([-500 500])
%if max(xph)>xmax; xlim([ | min(xph)<xmin; 
drawnow


toc;

%PLOT GRID CELLS DURING RETRIEVAL.
if plotGrids==1
figure('name','grid Retrieval','Numbertitle','off','Position',[250 50 numPhi*150 numF*150]); ax20=gca;

for ff=1:numF;
    f=baseF+multFF*ff;
    for pp=1:numPhi;
        for p2=1:numPhi;
plotSpikeX=xph(1,1:rrrT).*tspike3((ff-1)*numPhi^2+(pp-1)*numPhi+p2,1:rrrT);
plotSpikeY=yph(1,1:rrrT).*tspike3((ff-1)*numPhi^2+(pp-1)*numPhi+p2,1:rrrT);
axes(ax20);
subplot(numF,numPhi,(ff-1)*numPhi+pp,'replace'); ax20=gca;
plot(plotSpikeX,plotSpikeY,'r.');
xlim([xmin-1 xmax+1])
ylim([ymin-1 ymax+1])
drawnow
       end %End for p2=1:numPhi
    end %End for pp=1:numPhi
end %End for f=1:numF
end %End if plotGrids==1

%Plot Place Cells during retrieval.
if plotCells==1
figure('name','place cells Retrieval','Numbertitle','off','Position',[800 50 numP*100 numP*100]); ax21=gca;

for pp=1:numP^2;
plotPlaceX=xph(1,1:rrrT).*placeSpike3(pp,:);
plotPlaceY=yph(1,1:rrrT).*placeSpike3(pp,:);
axes(ax21);
subplot(numP,numP,pp,'replace'); ax21=gca;
plot(plotPlaceX,plotPlaceY,'r.');
xlim([xmin-1 xmax+1])
ylim([ymin-1 ymax+1])
end %End for pp
end %End if plotCells==1

if plotRet==1;
    
%PLOT PLACE VS TIME DURING RETRIEVAL.
axes(ax21b);
subplot(numPsizesTot*numPtestsTot,1,(numPsizes-1)*numPtestsTot+numPtests,'replace'); ax21b=gca;

numSING=numPlaceToPlot;
for sing=1:numSING;
    %for mm=1:length(timep(sing).spike);
    fpsPR=find(placeSpike3(sing,:));
    for mm=1:length(fpsPR);
        if colorSpikes==1;
            plot([fpsPR(mm) fpsPR(mm)],[sing sing-0.5],'Color',[rc(sing,:)],'LineWidth',1); hold on;
        elseif colorSpikes==0;
            plot([fpsPR(mm) fpsPR(mm)],[sing sing-0.5],'Color','k','LineWidth',1); hold on;
        end
        %plot([fps(mm) fps(mm)],[sing sing-0.5],'Color',[rc(sing,:)],'LineWidth',1); hold on;
        %plot([fps(mm) fps(mm)],[sing sing-0.5],'Color','k','LineWidth',1); hold on;
    end
end
if colorSpikes==1;
    set(gca,'Color','k');
elseif colorSpikes==0
    set(gca,'Color','w');
end
xlim([1 rrrT]);
ylim([0 numSING+0.5]);
drawnow

%PLOT HD CELLS VS SPACE during retrieval
%figure('name','HD cells Retrieval','Numbertitle','off','Position',[750 400 300 300]); ax19=gca;
axes(ax19);
subplot(numPsizesTot*numPtestsTot,1,(numPsizes-1)*numPtestsTot+numPtests,'replace'); ax19=gca; %This does one column of all.

plot(xph(1,2:rrrT),yph(1,2:rrrT),'-','LineWidth',1); %xlim([xmin+30 xmax+30]); ylim([ymin-60 ymax-60]); 
hold on

for m=2:20:rrrT-1;   %Plot the head direction cell vectors. %If want to plot those during encoding start m=1:
     denom=cos(h(1,1))*sin(h(2,1))-sin(h(1,1))*cos(h(2,1));
     hdx=((plothd3(1,m+1))*sin(h(2,1))-(plothd3(2,m+1))*sin(h(1,1))); %./(denom*wd*B);
     hdy=((plothd3(1,m+1))*-cos(h(2,1))+(plothd3(2,m+1))*cos(h(1,1))); %./(denom*wd*B);
    
      %hdx=hd3(1,m+1);
      %hdy=-cot(h(2,1))*hd3(1,m+1) + hd3(2,m+1)/sin(h(2,1));  %Translates back out of head direction coordinates.
      AA=[xph(1,m) xph(1,m)+hdx*0.5];  
      BB=[yph(1,m) yph(1,m)+hdy*0.5];   

      if sum(AA)~=0 && sum(BB)~=0;
          plot(xph(1,m), yph(1,m),'ko','LineWidth',1);
      plot(AA,BB,'k-','LineWidth',1); %AA,BB,'k-',
      xlim([xmin-10 xmax+10])
      ylim([ymin-10 ymax+10])
      hold on
      end
end
drawnow;

%PLOT HD ANGLE VERSUS TIME
axes(ax114);
tt=[1:rrrT];
plot(tt,plothd3(1,:));
xlim([1 rrrT]);
ylim([-20 20]);

%PLOT HD CELLS VS TIME DURING RETRIEVAL
axes(ax19b); %Plot HD retrieval versus time.
subplot(numPsizesTot*numPtestsTot,1,(numPsizes-1)*numPtestsTot+numPtests,'replace'); ax19b=gca;

fpsHDrrr().cell=[];
spikehd3=[];
spikehd3=zeros(length(h),length(hd3));
for ht=1:20:length(hd3); 
    maxhd3=mean(max(hd3)); %Create a spiking version of the HD activity in hd3.
    spikehd3(:,ht)=rand(length(h),1)>(maxhd3*ones(length(h),1)-hd3(:,ht))/maxhd3; 
    %Generate a spike when rand is greater than 1-hd3.
end
numSING=length(h);
for sing=1:numSING;
    fpsR=find(spikehd3(sing,:));
    fpsHDrrr(sing).cell=find(spikehd3(sing,:));
    for mm=1:length(fpsR);
%plot([fps(mm) fps(mm)],[sing sing-0.5],'Color','k','LineWidth',1); hold on;
        if colorSpikes==1;
            plot([fpsR(mm) fpsR(mm)],[sing sing-0.5],'Color',[rc(sing,:)],'LineWidth',1.8); hold on;
        elseif colorSpikes==0;
            plot([fpsR(mm) fpsR(mm)],[sing sing-0.5],'Color','k','LineWidth',1); hold on;
        end
    end
end
if colorSpikes==1;
    set(gca,'Color','k');
elseif colorSpikes==0
    set(gca,'Color','w');
end
xlim([1 rrrT]);
ylim([0 numSING+0.5]);
drawnow

end %End if plotRet==1

%PLOT distR
axes(ax22x);
subplot(numPsizesTot,numPtestsTot,(numPsizes-1)*numPtestsTot+numPtests,'replace'); ax22x=gca; 
    distR=sqrt((xph-xphE).^2 + (yph-yphE).^2);
    plot(1:rrrT,distR,'k','LineWidth',2); hold on; 
    ylim([0 100]);
    xlim([0 rrrT]);
    drawnow
meandistR=mean(distR);

exdistR(numPsizes,1)=numP^2; %arcnum; %numP;
exdistR(numPsizes,numPtests+1)=meandistR;


%PLOT X POS AND RUNNING SPEED FROM INVERSE TRANSFORM OF PHASE DURING
%RETRIEVAL
rspeed=zeros(1,rrrT);
for rt=3:rrrT;
    rspeed(1,rt)=sqrt((xph(rt)-xph(rt-1))^2 + (yph(rt)-yph(rt-1))^2)/dt;
end

if plotRet==1
axes(ax112);
tt=[1:rrrT];
plot(tt,yph(1,1:rrrT),tt,rspeed);
xlim([1 rrrT]);
ylim([ymin-10 ymax+10]);
end

    halfrun=0
    numCircRun(numPsizes,1)=numP;
    for xyt=1:rrrT-1
        if yph(xyt)<-5 && yph(xyt+1)>-5 && xph(xyt)<-10 && halfrun==0;
            numCircRun(numPsizes,numPtests+1)=numCircRun(numPsizes,numPtests+1)+0.5;
            halfrun=1;
        end
        if yph(xyt)>5 && yph(xyt+1)<5 && xph(xyt)>10 && halfrun==1;
            numCircRun(numPsizes,numPtests+1)=numCircRun(numPsizes,numPtests+1)+0.5;
            halfrun=0;
        end 
    end
    
    
    strengHD;
    numCircRun;
    %if datainput==5 && numCircRun(numPsizes,numPtests+1)>floor(strengHD*3)-1; stopWhile=10; end

    
        %end %End while numCircRun<2.0
    end %End numPsizes
end %End numPtests

%save WPGepBnAT WPG
%save -ascii WPGepBnAPascii WPG
%WPG
save EpBnAYhaftingEx exdistR
save -ascii EpBnAYhaftingEx exdistR

figure('name','mean dist ret-enc','Numbertitle','off','Position',[300 300 300 200]); ax22y=gca;
meanexdistR=mean(exdistR(:,2:length(exdistR(1,:)))');
plot(exdistR(:,1),meanexdistR,'LineWidth',2);
if arcnum>10
xlim([0 arcnum+10]);
elseif numP^2>10
    xlim([0 numP^2+10]);
end
ylim([0 100]);
    
fid=fopen('numCircRunAX','a');
fprintf(fid, '%g\t',numCircRun);
fprintf(fid, '\n');

fclose(fid);

toc

%END OF SCRIPT