% cpelaz.m
%
% format: cpelaz(fname,ID,DT,GLAT,GLNG,DECL)
%
% This program generates a summary of the GPS satellite
% positions in the sky and the associated scintillation 
% activity based off a computer S4 index.  The input
% data file is a .sum (summary) file.
%
% input:
%   fname - filename of the .sum file generated by sum.exe
%   ID - station ID, which is simply a character
%   DT - time difference in hours between UTC and LT (DT = UTC - LT)
%   GLAT - geographic latitude in degrees
%   GLON - geographic longitude in degrees
%   DECL - magnetic declination in degrees (note: west is negative)
%
% example parameters:
%
% fname = '010118X0.sum';
% ID = 'X';
% DT = 3;
% GLAT = 42.44;   
% GLON = -76.48;
% DECL = -21;
%
%
%
% original:  h. kil
% update:    b. ledvina  9.17.01  v2.0


function cpelaz2(fname,ID,DT,GLAT,GLNG,DECL)

%***************************************************************
%fname='010118X0.sum'; 
%DT=3;        %time difference between UT and LT: DT=UT-LT 
%ID='X';      %station ID
%DECL=-21;    %magnetic declination (east is positive)
%GLNG=-45.0;  %geographic longitude
%GLAT=-22.4;  %geographic latitude
%MLAT=-13;    %magnetic latitude


      %*** Plotting parameters ***
       xl=0.2;    % the following four paremeters
       yb=0.3;    %   are used to indicate
       xr=8;      %   paper positions
       yr=10;

       tt=linspace(0,2*pi);
       xx1=90*cos(tt); yy1=90*sin(tt);
       xx2=60*cos(tt); yy2=60*sin(tt);
       xx3=30*cos(tt); yy3=30*sin(tt);
       xy=90/sqrt(2);
       decx=90/sqrt(tan(pi/2-DECL*pi/180)^2+1);
       decy=tan(pi/2-DECL*pi/180)*decx;
       
       SUBPOS=[0.08 0.73 0.22*yr/xr 0.22;0.08 0.50 0.22*yr/xr 0.22;
               0.08 0.27 0.22*yr/xr 0.22;0.08 0.04 0.22*yr/xr 0.22;
               0.36 0.73 0.22*yr/xr 0.22;0.36 0.50 0.22*yr/xr 0.22;
               0.36 0.27 0.22*yr/xr 0.22;0.36 0.04 0.22*yr/xr 0.22;
               0.64 0.73 0.22*yr/xr 0.22;0.64 0.50 0.22*yr/xr 0.22;
               0.64 0.27 0.22*yr/xr 0.22;0.64 0.04 0.22*yr/xr 0.22]; 
            
      LTLAB='18-19-20-21-22-23-00-01-02-03-04-05-06';
      %**********************************************************
            
  
%************************************************************      
fid=fopen(fname,'r');  %open summary a file
    
    %*** READ HEADER ***
      LOCA=fgets(fid);
      if length(LOCA) >= 30
        CITY = LOCA(1:30);
      else
        CITY = LOCA(1:end);
      end
      COMM=fgets(fid);
      dum1=find(LOCA==13 | LOCA==10);
        if (~isempty(dum1))
           LOCA=LOCA(1:(min(dum1)-1));
        end
      dum2=find(COMM==13 | COMM==10);
        if (~isempty(dum2))
           COMM=COMM(1:(min(dum2)-1));
        end

      HD=fscanf(fid,'%f %f %f %d %d %d %d %d %d',[9]); 
         Xo=HD(1); Yo=HD(2); Zo=HD(3); YE  =HD(4); MO=HD(5);
         DA=HD(6); HO=HD(7); MI=HD(8); Nrec=HD(9);
         
     %*** READ DATA ***   
      i=ones(32,1); UT=zeros(32,800); 
      Xsat=UT; Ysat=UT; Zsat=UT; S4=UT; EL=UT; AZ=UT; WBP=ones(32,800);   
     
      stop=0;
      OLDTIME=9999;
      HHMM=1900+DT*100;

      UUU=1;
      
    while(stop==feof(fid))

       TIME=fscanf(fid,'%f %d %f',[3]);
       if (isempty(TIME)==1), break,end

       if(TIME(1) == OLDTIME | TIME(1)<1800+DT*100)
          for Nsat=1:TIME(2)
            ENTRY=fscanf(fid,'%d %f %f %f %d %d %d %f %d %d',[10]);
          end
       else  
          for Nsat=1:TIME(2)
            ENTRY=fscanf(fid,'%d %f %f %f %d %d %d %f %d %d',[10]);
              if (ENTRY(8)<2 & ENTRY(8)>0)
                PRN                = ENTRY(1); 
                Xsat(PRN,i(PRN,1)) = ENTRY(2);
                Ysat(PRN,i(PRN,1)) = ENTRY(3);
                Zsat(PRN,i(PRN,1)) = ENTRY(4);  
                S4  (PRN,i(PRN,1)) = ENTRY(8);
                UT  (PRN,i(PRN,1)) = TIME (1);
                
                WBP (PRN,i(PRN,1)) = 10*log10(ENTRY(6));

		elaz=sum_elaz(Xo,Yo,Zo,ENTRY(2),ENTRY(3),ENTRY(4));

                EL(PRN,i(PRN,1))   = elaz(1);
                AZ(PRN,i(PRN,1))   = elaz(2);
                       i(PRN,1)    = i(PRN,1)+1;
                
		end      
          end
       end
      OLDTIME = TIME(1);  
       
      %*** PLOT ****
       if (TIME(1)>=HHMM)
  
          %*** MAKE CIRCLES ***
           subplot('position',SUBPOS(UUU,:))
           plot(xx1,yy1),axis([-100 100 -100 100])
            hold on
           plot(xx2,yy2)
            hold on
           plot(xx3,yy3)
            hold on
           plot([-90 90],[0 0],':')
            hold on
           plot([0 0],[-90 90],':')
            hold on
           plot([-xy xy],[-xy xy],':')
            hold on
           plot([-xy xy],[xy -xy],':')
            hold on
           plot([-decx decx],[-decy decy],'--')
              arrow(decx,-decx,decy,-decy,2,7,15);
            hold on
              arrow(-decx,decx,-decy,decy,2,7,15);
            hold on
         
           %*** SATELLITE PATH ***
                   
             for PRN=1:32     
                  NI=i(PRN,1)-1;
               if (NI>30) 
                     x=ones(NI,1);
                     y=ones(NI,1);
                     
                 %*** Convert elevation and azimuth to x and y coordinates. 
                  for jj=1:NI
                      eldeg=EL(PRN,jj);
                      azrad=AZ(PRN,jj)*pi/180;
                      x(jj)=(90-eldeg)/sqrt(1+cot(azrad)^2);
                    if (azrad>pi)
                      x(jj)=-x(jj);
                    end
                    % changing code for when azrad = 0;
                    % brent m ledvina
                    %if(azrad == 0)
                     %  y(jj) = x(jj);
                    %else
                      y(jj)=x(jj)*cot(azrad);
                    %end 
                     
                  end 
                      plot(x,y) 
                       hold on
                      arrow(x(NI),x(NI-10),y(NI),y(NI-10),1,6,5);    

                 %*** make the satellite path width proportion to S4 index. 
                  dx=diff(x);
                  dy=diff(y);
                  
                  jjj=find(dx~=0);
                  sjjj=size(jjj);
                  
                      xys4=[];
                 if (sjjj(1)>30)    
                   for nn=1:3     % 3-point mean
                      for mm=2:NI-1 
                         S4(PRN,mm)=(S4(PRN,mm-1)+S4(PRN,mm)+S4(PRN,mm+1))/3;
                      end   
                   end
                         ss=find(S4(PRN,1:NI)<0.1);
                         sizess=size(ss);
                         if (sizess(1)>0)
                            S4(PRN,ss)=0;
                         end   
                   for kk=1:sjjj(1)
                      theta=atan(dy(jjj(kk))/dx(jjj(kk)));
                      R=[cos(theta) -sin(theta);sin(theta) cos(theta)];
            
                      delta1=R*[0;  S4(PRN,jjj(kk))*10];
                      delta2=R*[0; -S4(PRN,jjj(kk))*10];
 
                      xu=x(jjj(kk))+delta1(1);
                      yu=y(jjj(kk))+delta1(2);
                      xd=x(jjj(kk))+delta2(1);
                      yd=y(jjj(kk))+delta2(2);
             
                      xys4=[xys4;xu yu xd yd];
                    end
                       hold on
                      plot(xys4(:,1),xys4(:,2))
                       hold on
                      plot(xys4(:,3),xys4(:,4))
                       hold on
                        if (dx(sjjj(1))<0)
                           theta=theta+pi;
                        end
                        
                           slp  = (y(2)-y(1))/(x(2)-x(1));
                           PRNx = -(x(2)-x(1))*10+x(1);
                           PRNy = (PRNx-x(1))*slp+y(1);
                           if (abs(PRNx-x(1))>15 | abs(PRNy-y(1))>15)
                              PRNx=x(1); PRNy=y(1);
                           end   
                        text(PRNx,PRNy,num2str(PRN),'HorizontalAlignment','center',...
                           'Fontsize',8)
       
                 end   % end plotting satellite path     
               end   %end  if (NI>10)      
             end  % end for PRN=1:NI           

             hold off 
             axis off                

              text(30,90,LTLAB(3*UUU-2:3*UUU+2),'Fontsize',11,...
                    'VerticalAlignment','bottom','HorizontalAlignment','left'  )
           if (UUU==1)
              text(-90,138,[CITY,' (',ID,')'],'Fontsize',15)
              text(390,138,[num2str(DA),'.',num2str(MO),'.',num2str(YE)],...
                                'HorizontalAlignment','Left','Fontsize',16);

           end   
           if (UUU==1 | UUU==5 | UUU==9)
              text(0,92,'\fontsize{10}0\circ','VerticalAlignment','Bottom',...
                                               'HorizontalAlignment','Center' )
           end  
           if (UUU>=9 & UUU<=12)
              text(93,0,'\fontsize{10}90\circ','VerticalAlignment','middle',...
                                               'HorizontalAlignment','left' )
           end
           if (UUU==4 | UUU==8 | UUU==12)
              text(0,-92,'\fontsize{10}180\circ','VerticalAlignment','Top',...
                                               'HorizontalAlignment','Center' )
           end
           if (UUU<=4)
              text(-93,0,'\fontsize{10}270\circ','VerticalAlignment','middle',...
                                               'HorizontalAlignment','right' )
           end
    
       %********************************************************************

                
          HHMM=TIME(1)+100;
          UUU=UUU+1;
          i=ones(32,1);   
     
       end  %end if (HHMM >)    
       
        if (UUU>12),break,end
 end  %end while 
        set(gcf,'paperposition',[xl yb xr yr]);     
        fclose(fid);






































































































