%This function outputs a simulated survey (Output.X with columns [site number,
%time, x-coordinate, y-coordinate]) and the resulting variance-covariance matrix
%with observation error (Output.VarCov).
%Function arguments: TotalSamples = total number of samples (ST), LandSize = vector of length 2 representing
%the dimensions of the landscape, SurveyYears = number of years in which surveys take place, Years = length
%of monitoring period, Sigma = standard deviation of environmental variation (sigma), OneMinusGamma = 1 - gamma, 
%Rho = rho, SigmaObs = standard deviation of observation error (sigma_obs).
%Example: getsurveyo(100,[10 10],2,10,0.5,0.5,0.5,0.5);

function Output = getsurveyo(TotalSamples,LandSize,SurveyYears,Years,Sigma,OneMinusGamma,Rho,SigmaObs)

%check that Years >= SurveyYears
if (Years < SurveyYears)
    error('years surveyed is greater than the number of years');
end

%check that TotalSamples is divisible by SurveyYears 
if (mod(TotalSamples,SurveyYears) ~= 0)
    error('total sample size must be divisible by the number of survey years');
end

%set up output matrics
X = ones(TotalSamples,4);
Cov = zeros(TotalSamples,TotalSamples);

%set up matrix to hold site number values
SiteNumbers = zeros(TotalSamples,1);

%get the current sample number
SampleNumber = 1;

%get the number of sites
Sites = TotalSamples / SurveyYears;

%loop through the sites and allocate values to X
for (Site = 1:Sites)
   %reset within site number and sites left within site
   WithinSiteNum = 0;
   WithinSiteLeft = SurveyYears;
   YearsLeft = Years;
   
   %get random lat and long for sites 
   Lat = rand(1) * LandSize(1);
   Long = rand(1) * LandSize(2);
   
   %loop through surveys and record value for X and Site Numbers
   for (Sample = SampleNumber:(SampleNumber + SurveyYears - 1))
       SiteNumbers(Sample,1) = Site;
       %get the current time for this survey 
       if (Sample == SampleNumber)
            Time = 1;
            YearsLeft = Years - Time;
       elseif (Sample == (SampleNumber + SurveyYears - 1))
            Time = Years;
            YearsLeft = 0;
       else
            Time = Time + round((YearsLeft / WithinSiteLeft));
            YearsLeft = Years - Time;
       end     
       
       %record the time and location
       X(Sample,1) = Sample;
       X(Sample,2) = Time;
       X(Sample,3) = Lat;
       X(Sample,4) = Long;
             
       %increment within site number
       WithinSiteNum = WithinSiteNum + 1;
       WithinSiteLeft = WithinSiteLeft - 1;
   end
      
   %increment sample number
   SampleNumber = SampleNumber + SurveyYears;
end

%loop through samples and populate the variance-covariance matrix
for (Sample1 = 1:TotalSamples)
   for(Sample2 = 1:TotalSamples)
      if (SiteNumbers(Sample1,1) == SiteNumbers(Sample2,1))
          if (X(Sample1,2) == X(Sample2,2))
              Cov(Sample1,Sample2) = ((Sigma ^ 2) / (1 - (OneMinusGamma ^ 2))) + (SigmaObs ^ 2);
          else
              Cov(Sample1,Sample2) = ((Sigma ^ 2) / (1 - (OneMinusGamma ^ 2))) * (OneMinusGamma ^ abs(X(Sample1,2) - X(Sample2,2)));
          end
      else
          if (X(Sample1,2) == X(Sample2,2))
              Cov(Sample1,Sample2) = ((Sigma ^ 2) / (1 - (OneMinusGamma ^ 2))) * (Rho ^ sqrt(((X(Sample1,3) - X(Sample2,3)) ^ 2) + ((X(Sample1,4) - X(Sample2,4)) ^ 2)));
          else
              Cov(Sample1,Sample2) = ((Sigma ^ 2) / (1 - (OneMinusGamma ^ 2))) * (Rho ^ sqrt(((X(Sample1,3) - X(Sample2,3)) ^ 2) + ((X(Sample1,4) - X(Sample2,4)) ^ 2))) * (OneMinusGamma ^ abs(X(Sample1,2) - X(Sample2,2)));
          end
      end
   end
end

Output.X = X;
Output.VarCov = Cov;