%returns average orbit mesured from devices devdescr on basis of given number (navg) of measuremns
function [hmean,vmean,currmean] = genorbit()

 global spos thpat testmode ncalls scalings

 l = length(spos);


 fprintf(1,'ncalls = %d\n',ncalls);
  
 errfact = 10;
 if (ncalls == 0)
 
    fprintf(1,'Generating the reference orbits of length %d\n',l);

    hmean(:,1) =  normrnd(0,3,[l 1]);
    vmean(:,1) =  normrnd(0,1 ,[l 1]);

    hmean(:,2) =  poissrnd(100,[l 1])/(100 *errfact);
    vmean(:,2) =  poissrnd(100,[l 1])/(1000*errfact);

    currmean = zeros(l,2);
    currmean(:,1) = -4;
    currmean(:,2) = 0.2;

    ncalls = ncalls + 1;
    return;
 
 else  
 
 
    global referenceH referenceV referenceC 

    scaling = 2.5;
    hmean(:,1)    = thpat(:,1)*scalings(ncalls)*scaling;% generates the pattern prop to the the patterb
    vmean(:,1)    = zeros(l,1);
    currmean(:,1) = zeros(l,1);
    currmean(:,1) = -4;

    hmean(:,1)  = hmean(:,1) + normrnd(0,3.0,[l 1])/errfact; 
    vmean(:,1)  = vmean(:,1) + normrnd(0,0.1,[l 1])/errfact; 


    hmean(:,1)  = hmean(:,1) + referenceH(:,1); 
    vmean(:,1)  = vmean(:,1) + referenceV(:,1); 


    hmean(:,2)    = poissrnd(100,[l 1])/(100 *errfact); % H error bars
    vmean(:,2)    = poissrnd(100,[l 1])/(1000*errfact);% V error bars
    currmean(:,2) = poissrnd(50,[l 1])/ (100*errfact);  % C error bars

    ncalls = ncalls + 1;                                                

    return;                                                             


 end                                                                             
                                                                                






