function part2sim2powersin(K,L,N,sigma,numobs)
  %L is here assumed known
  angl = pi/4;
  d = 1;
  lambda = 10;
  
  I2 = lambda*log( (1+sin(angl))/(1-sin(angl)) ) / ( 4*d*(angl^2) );
  I3 = (lambda^2)*tan(angl)/(4*(d^2)*(angl^3));
  
  P = zeros(L,L);
  step = floor(L/3);
  for (k=1:step)
    P(k,k) = 0.5;
    P(k+step,k+step) = 1;
    P(k+2*step,k+2*step) = 1.5;
  end
  
  c1 = N/K;
  c2 = L/N;
  c3 = L/K;
  eigs1 = []; eigs2 = []; eigs3 = [];
  obsvs = linspace(1,numobs,100);
  obsvs = floor(obsvs);
  for (obsnr=obsvs)
    obsnr
    samplecovsum1 = 0;
    samplecovsum2 = 0;
    samplecovsum3 = 0;
    for (thisobs=1:obsnr)
      V = zeros(N,L);
      for (k=1:L) 
        theta = 2*angl*(rand(1)-0.5);
        V(:,k) = ((exp(2*pi*j*sin(theta)*d/lambda)).^(0:(N-1)))';
      end
      V = (1/sqrt(N))*V;
      alpha = randn(L,K);
      nmatrix = sigma*randn(N,K);
      Y = V*sqrt(P)*alpha + nmatrix;
      W = (1/K)*Y*Y';%sample covariance matrix sample
      samplecovsum1 = samplecovsum1 + (1/N)*trace(W);
      samplecovsum2 = samplecovsum2 + (1/N)*trace(W^2);
      samplecovsum3 = samplecovsum3 + (1/N)*trace(W^3);
    end
    W1 = samplecovsum1/obsnr;
    W2 = samplecovsum2/obsnr;
    W3 = samplecovsum3/obsnr;
                 
    P1 = (W1 -sigma^2)/c2;
    P2 = ( W2     -      ((c2^2)*I2 + c2*c3)*(P1)^2 ...     
                  -      2*(sigma^2)*(c2+c3)*P1 ...    
                  -      (sigma^4)*(1+c1) ...
         )/c2;
    P3 = (W3      -      (3*(c2^2)*I2 + 3*c2*c3)*P1*P2 ...
                  -      ( (c2^3)*I3 + 3*(c2^2)*c3*I2 + c2*(c3^2) )*(P1^3) ...    
                  -      3*(sigma^2)*( (1+c1)*c2 )*P2 ...    
                  -      3*(sigma^2)*( (1+c1)*(c2^2)*I2 + c3*(c3+2*c2) )*(P1^2) ...    
                  -      3*(sigma^4)*( c1^2+3*c1+1)*c2*P1 ...     
                  -      (sigma^6)*(c1^2+3*c1+1) ...      
         )/c2;
    [P1 P2 P3]
    eigs = abs(eigenvalsfrommoments([P1 P2 P3]));
    eigs = sort(eigs);
    eigs1 = [eigs1 eigs(1)];
    eigs2 = [eigs2 eigs(2)];
    eigs3 = [eigs3 eigs(3)];
  end  
  
  plot(obsvs,eigs1,'k+',obsvs,eigs2,'kx',obsvs,eigs3,'ko');
  xlabel('Number of observations');
  legend('First eigenvalue','Second eigenvalue','Third eigenvalue');