global tumin mu radiusearthkm xke j2 j3 j4 j3oj2  
format longG;
% mode
popis=true;
kde='n';
rad=180/pi;
vybrano=false;
twopi=2*pi;
outfile2= fopen('uloha2.a.log', 'wt');
GEAkonstanty;
constastro;
c1=system('wget -q -O geodetic.txt http://www.celestrak.com/NORAD/elements/geodetic.txt');
infile = fopen('geodetic.txt', 'r');
if (infile == -1)
       fprintf(1,'Failed to open file: %s\n', 'geodetic.txt');
       return;
end
% fi=input('zadej sirku. (deg) ')/rad;
% lambda=input('zadej delku (E) (deg) ')/rad;
fi=50.1/rad;
lambda=14.25/rad;
beta=10/rad;
obs=sf2xy(fi,lambda,'r2')*RZ;
krok=60;
krokden=krok/86400;
while (~feof(infile))
   jmenop=fgetl(infile);
   vybrano=input(jmenop)
   if(!vybrano)
      fgetl(infile);
      fgetl(infile);
      continue
   end   
   longstr1 = fgetl(infile);
   while ( (longstr1(1) == '#') && (feof(infile) == 0) )
        longstr1 = fgetl(infile);
   end
   if (feof(infile) == 0)
        
       longstr2 = fgetl(infile);
       i=1;  
       if(vybrano)
            [satrec] = readtwoline(longstr1, longstr2);
            jmenoc=satrec.intldesg;
            year=satrec.epochyr+2000;
            mjd0=gr2jd(1,1,year,'true')-1;
            mjdep=mjd0+satrec.epochdays;
            n0=satrec.no/86400*2*pi;
            Perioda=86400/satrec.no;
            elementy(1)=nthroot(GMZ/n0^2,3);
            elementy(2)=satrec.ecco;
            elementy(3)=satrec.inclo;
            elementy(4)=satrec.nodeo;
            elementy(5)=satrec.argpo;
            elementy(6)=n0;
            M00=satrec.mo;
            elementy(7)=M00;
            fprintf(outfile2,'# %s %10.2f %8.6f %7.2f %7.2f %7.2f %7.2f\n',jmenop,elementy(1:5),M00*rad);
         end
         if(vybrano)
           break
         end
      end
 end
  if(vybrano)
%  urceni viditelnosti druzice beta stupnu nad obzorem, ze 
%  zadane hlavni poloosy drahy
    Rla=RZ/elementy(1);
    tbeta=tan(beta);
%  2 reseni kvadraticke rovnice     
    sina21=(-Rla*tbeta+sqrt(1-Rla^2+tbeta^2))/(1+tbeta^2);
    sina22=(-Rla*tbeta-sqrt(1-Rla^2+tbeta^2))/(1+tbeta^2);
    al2a=asin(sina21);
    al2b=asin(sina22);
% vyber realneho reseni
    if((al2a<1.4) &  (al2a > 0)) % geostat je maximum 81 stupnu
      alfa2=al2a;
    else
      alfa2=al2b;
    end
    dvid=RZ*alfa2; % vzdálenost, kdy je druzie jeste nad obzorem.
 % delta lambda - pootoceni Země za jeden oběh družice
    dlambda=OmegaZ*Perioda; %(dlambda v rad)
    disp(['dlambda=',num2str(dlambda*rad)]);
   
    disp('cas pro ktery pocitame');
    disp(['rok=',num2str(year)]);
    mesic=input('mesic=?')
    den = input('den=?')
    hodina=input('hodina=?')
    minuta=input('minuta=?')
    uth=(hodina+minuta/60);
    utd=uth/24;
    mjd=gr2jd(den,mesic,year,true)+utd;
    obzsour=[2,0,0];
    ii=0;
    dtt=1./1440.;
    while (obzsour(1)*rad >= 1 )
     dth=ii/60
     mjdx=mjd+ii*dtt;
     utc=uth+dth
     obzsour=topozelem20(mjdx,mjdep,elementy,obs)
     fprintf(outfile2,'%2d%3d%5d%8.4f',den,mesic,year,utc);
     fprintf(outfile2,'%9.4f%9.4f%12.1f\n',obzsour(1)*rad,obzsour(2)*rad,obzsour(3));
     disp(num2str(mjdx)) 
     disp(num2str(obzsour(1)*rad))
     disp(num2str(obzsour(2)*rad));
     disp(num2str(obzsour(3)));
     ii = ii + 1;
   end 
   fclose(outfile2);
end
