% 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);