function Output=GravitySpec(mapO,thRateRatio,PhasePlot,Fplot)

%%%%%%mapO: Spectrum map (Fre*Phase)
%%%%%%Gravity Threshold of Place Field was selected thRateRatio*PeakFiringRate
%%%%%%Place field was defined as >thNum of pixel with firing rate>thRateRatio*PeakFiringRate
PeakPhase=[];
PeakFre=[];
GravityPhase=[];
GravityFre=[];

if length(size(mapO))==2
    [PeakPhase,PeakFre,GravityPhase,GravityFre,GravNum,GravityPhaseStd,GravityFreStd]=GravitySpecMat(mapO,thRateRatio,PhasePlot,Fplot);   
elseif length(size(mapO))==3
    for j=1:size(mapO,3)
    [PeakPhase(j),PeakFre(j),GravityPhase(j),GravityFre(j),GravNum(j),GravityPhaseStd(j),GravityFreStd(j)]=GravitySpecMat(squeeze(mapO(:,:,j)),thRateRatio,PhasePlot,Fplot);   

    end
    elseif length(size(mapO))==4
    for j=1:size(mapO,3)
        for i=1:size(mapO,4)
            [PeakPhase(j,i),PeakFre(j,i),GravityPhase(j,i),GravityFre(j,i),GravNum(j,i),GravityPhaseStd(j,i),GravityFreStd(j,i)]=GravitySpecMat(squeeze(mapO(:,:,j,i)),thRateRatio,PhasePlot,Fplot);   
        end
    end


else

end
Output.PeakPhase=PeakPhase;
Output.PeakFre=PeakFre;
Output.GravityPhase=GravityPhase;
Output.GravityFre=GravityFre;
Output.GravNum=GravNum;
Output.GravityPhaseStd=GravityPhaseStd;
Output.GravityFreStd=GravityFreStd;



function [PeakPhase1,PeakFre1,GravityPhase1,GravityFre1,GravNum,GravityPhaseStd1,GravityFreStd1]=GravitySpecMat(map1,thRateRatio1,PhasePlot1,Fplot1)

   thNum=30;
   Cell=PlaceFieldFind2D(map1,thRateRatio1,thNum);
   GravX1=[];
   GravY1=[];
   for i=1:length(Cell)
       GravNum(i)=length(Cell(i).IndX); 
   end
   [GravNum,II1]=max(GravNum);
       GravX1=Cell(II1).IndX(:);
       GravY1=Cell(II1).IndY(:);

   
   PeakRate1=max(max(map1));
   [GravX1,GravY1,w1]=find(map1>=PeakRate1*thRateRatio1);
%    Weighted=map1(GravX1);
%    Weighted=Weighted-min(Weighted);
%    Weighted=Weighted/sum(Weighted);
   
   for i=1:length(GravX1)
       Weighted(i)=map1(GravX1(i),GravY1(i));
   end
   Weighted=Weighted-min(Weighted);
   Weighted=Weighted/sum(Weighted);

   GravNum=length(GravX1);
   [PeakX1,PeakY1]=find(map1>=PeakRate1);
   PeakPhase1=PhasePlot1(PeakX1);
   PeakFre1=Fplot1(PeakY1);
   
   kkk=Fplot1(GravY1);
   GravityFre1=sum(kkk(:).*Weighted(:));

   temp=PhasePlot1(GravX1);
   m=zeros(length(PhasePlot1),1);
   for i=1:length(PhasePlot1)
       kk=find(GravX1==i);
       if ~isempty(kk)
       m(i)=sum(Weighted(kk));
       end
   end
%    m=m.*Weighted;
%    GravityPhase1=circ_mean(temp(:));
    GravityPhase1 = circ_mean(PhasePlot1(:), m(:));
    GravityPhaseStd1 = circ_std(PhasePlot1(:), m(:));
    
    PopulationNum=100;
    SampleK=PopulationNum*Weighted(:);
    GravityFreStd1=sqrt(sum(SampleK(:).*(kkk(:)-GravityFre1).^2)/PopulationNum);
    


