No BSD License  

Highlights from
Skyplot1001

image thumbnail
from Skyplot1001 by Markus Penzkofer
Virtual Planetarium showing stars, the planets, the sun and the moon.

pplotsky(arg)
function pplotsky(arg)

% function pplotsky(arg)
%
% Virtual Planetarium with Graphical User Interface
% - program
%
% 05.05.04 M.Penzkofer

% global variables: Handles
global H_LAMBDA H_PHI H_ZONE H_UT H_NZ
global H_YEAR H_MONTH H_DAY H_HOUR H_MJD
global H_CHB H_PICK
global NUM
global anz_fk5 anz_ppm num recall decall mag rec dec az h
global T recP decP recM decM azP hP azM hM

% location of observation:
% geogr. longitude (eastern neg.)
lambda = str2num(get(H_LAMBDA,'String'));
% geogr. latitude
phi = str2num(get(H_PHI,'String'));
% time of zone - UT (in h)
zone = str2num(get(H_ZONE,'String'));

% date / time :
% year, month, day, hour
jahr = str2num(get(H_YEAR,'String'));
monat = str2num(get(H_MONTH,'String'));
tag = str2num(get(H_DAY,'String'));
stunde = str2num(get(H_HOUR,'String'));

% show names of constellations:
showconst = get(H_CHB,'Value');

% ---------------------------------------------------------------------------------------------------------

% check of the input

% completeness
if (isempty(lambda) | isempty(phi) | isempty(zone) | isempty(jahr) | isempty(monat) | isempty(tag) | isempty(stunde))
   ErrorString = 'Enter all input values!';
   disp(ErrorString);
   ErrMsg(ErrorString);
   return
end

% format
if (frac(jahr) ~= 0 | frac(monat) ~= 0 | frac(tag) ~= 0)
   ErrorString = 'Enter integer values for the date!';
   disp(ErrorString);
   ErrMsg(ErrorString);
   return
end

% location of observation
if (abs(lambda) > 180)
   ErrorString = 'Enter LONGITUDE between -180 and +180 Degrees!';
   disp(ErrorString);
   ErrMsg(ErrorString);
   set(H_LAMBDA,'String','0');
   return
end
if (abs(phi) > 90)
   ErrorString = 'Enter LATITUDE between -90 and +90 Degrees!';
   disp(ErrorString);
   ErrMsg(ErrorString);
   set(H_PHI,'String','0');
   return
end
if (zone < -12 | zone > 12)
   ErrorString = 'Enter UT difference between -12 and +12 h!';
   disp(ErrorString);
   ErrMsg(ErrorString);
   set(H_ZONE,'String','0');
   return
end

% date / time
if (jahr < -4713 | jahr > 10000)
   ErrorString = 'Enter YEAR between -4713 and +10000!';
   disp(ErrorString);
   ErrMsg(ErrorString);
   set(H_YEAR,'String','0');
   return
end
if (monat < 1 | monat > 12)
   ErrorString = 'Enter MONTH between 1 and 12!';
   disp(ErrorString);
   ErrMsg(ErrorString);
   set(H_MONTH,'String','1');
   return
end
if (tag < 1 | tag > 31)
   ErrorString = 'Enter DAY between 1 and 31!';
   disp(ErrorString);
   ErrMsg(ErrorString);
   set(H_DAY,'String','1');
   return
end
if (stunde < 0 | stunde > 23.99)
   ErrorString = 'Enter HOUR between 0.00 and 23.99!';
   disp(ErrorString);
   ErrMsg(ErrorString);
   set(H_HOUR,'String','0');
   return
end


% ---------------------------------------------------------------------------------------------------------

% input files:
ffk5 = 'starsfk5.dat';
fppm = 'starsppm.dat';
fpic = 'constellations.dat';

% read input files
eval(['load ' ffk5 ' -ascii']);
eval(['load ' fppm ' -ascii']);

% constellations
fid=fopen(fpic);
	c = 0;
	name = fscanf(fid,'%s',[1 1]);
	while ( isequal(name,'') == 0 )
      c = c+1;
      indus = find(name == '_');
      if (indus > 0)
         name(indus) = ' ';
      end
      sternbilder(c).name = name;
   	recC(c) = fscanf(fid,'%f',[1 1])*15.0;
      decC(c) = fscanf(fid,'%d',[1 1]);
      name = fscanf(fid,'%s',[1 1]);
   end
   anz_bilder = c;
fclose('all');

% starts: FK5/PPM number, rectascension, declination, proper motion and magnitude, J2000
sterne1 = starsfk5;
sterne2 = starsppm;
num = [sterne1(:,1); sterne2(:,1)];
recall(:,1:4) = [sterne1(:,2:5); sterne2(:,2:5)];
decall(:,1:4) = [sterne1(:,6:9); sterne2(:,6:9)];
mag = [sterne1(:,10); sterne2(:,10)];

% number of stars
anz_fk5 = size(sterne1,1);
anz_ppm = size(sterne2,1);
anz_sterne = anz_fk5+anz_ppm;

% time zone
zone_num = -round(lambda/15);
set(H_NZ,'String',num2str(zone_num));
% UT
UT = stunde-zone;
set(H_UT,'String',num2str(UT));
% Julian date
mjdat = mjuldat(jahr,monat,tag,stunde)-zone/24.0;
set(H_MJD,'String',num2str(mjdat+2400000.5));
% Julian century in respect of J2000
T = jjh2000(mjdat,1);
% equinox of catalogue
TEQX = 0.0;

% rectascension, declination in decimal degree
for i = 1:anz_sterne
   rec(i) = (gms2alt(recall(i,1),recall(i,2),recall(i,3))+gms2alt(0,0,recall(i,4)*(T-TEQX)))*15.0;
   dec(i) = gms2alt(decall(i,1),decall(i,2),decall(i,3))+gms2alt(0,0,decall(i,4)*(T-TEQX));
end

% matrix for precession and nutation
% from mean to true equinox
A = pnmatrix(TEQX,T);
% aberration
[vx,vy,vz] = aberrat(T);

% constellations
for c = 1:anz_bilder
	% cartesian coordinates
   [xC(c),yC(c),zC(c)] = pol2cart(1,decC(c),recC(c));
   % precession
   [xC(c),yC(c),zC(c)] = precart(A,xC(c),yC(c),zC(c));
	% aberration
   xC(c) = xC(c)+vx; yC(c) = yC(c)+vy; zC(c) = zC(c)+vz;
	% new polar apparent coordinates
	[rC(c),decC(c),recC(c)] = cart2pol(xC(c),yC(c),zC(c));
end
% stars
for i = 1:anz_sterne
	% cartesian coordinates
   [x(i),y(i),z(i)] = pol2cart(1,dec(i),rec(i));
   % precession
   [x(i),y(i),z(i)] = precart(A,x(i),y(i),z(i));
	% aberration
   x(i) = x(i)+vx; y(i) = y(i)+vy; z(i) = z(i)+vz;
	% new polar apparent coordinates
	[r(i),dec(i),rec(i)] = cart2pol(x(i),y(i),z(i));
end

% planets incl. sun
planetenStr = str2mat(' Mercury ',' Venus ',' SUN ',' Mars ',' Jove ',' Saturn ',' Uranus ',' Neptune ',' Pluto ');
% geozentr. equator. positions incl. sun
[recP,decP] = planequ(T);

% moon
[recM,decM,rM] = moonequ(T);
 
% local mean sidereal time
LMST = lmst(mjdat,lambda);
% hour angle for stars
tau = (15*LMST-rec);
% horizon system stars
for i = 1:anz_sterne
	[h(i),az(i)] = equ2hor(dec(i),tau(i),phi);
end

% hour angle for planets
tauP = (15*LMST-recP);
% horizon system planets
for p = 1:9
	[hP(p),azP(p)] = equ2hor(decP(p),tauP(p),phi);   
end

% hour angle for moon
tauM = (15*LMST-recM);
% horizon system mond
[hM,azM] = equ2hor(decM,tauM,phi);

% hour angle for constellations
tauC = (15*LMST-recC);
% horizon system constellations
for c =1:anz_bilder
	[hC(c),azC(c)] = equ2hor(decC(c),tauC(c),phi);
end

% selection of stars by magnitude
ind00 = find(-1.5 < mag & mag <= 0.0);
ind01 = find(-0.0 < mag & mag <= 1.5);
ind02 = find(1.5 < mag & mag <= 3.0);
ind03 = find(3.0 < mag);

% background dependent on elevation of sun
if (hP(3) < -18)
	backCol = [0 0 0];
elseif (-18 <= hP(3) & hP(3) <= 0)
   backCol = [0.5 0.5 0.5];
elseif (0 < hP(3))
   backCol = [0 0 1];
end

% flag done
NUM = NUM+1;
% clear axis
if (NUM > 1)
   cla
end
% axes and labeling
if (isequal(arg,'P') == 1)
   % panorama
   set(gca,'Position',[0.05 0.25 0.90 0.30],'Color',backCol);
   axis([0 360 0 70]);
   HT1=text(0,-10,'S','FontSize',14,'HorizontalAlignment','Center');
	HT2=text(90,-10,'W','FontSize',14,'HorizontalAlignment','Center');
	HT3=text(180,-10,'N','FontSize',14,'HorizontalAlignment','Center');
	HT4=text(270,-10,'E','FontSize',14,'HorizontalAlignment','Center');
	HT0=text(360,-10,'S','FontSize',14,'HorizontalAlignment','Center');
elseif (isequal(arg,'W') == 1)
   % West
   set(gca,'Position',[0.05 0.07 0.90 0.55],'Color',backCol);
   axis([0 180 0 70]);
	HT1=text(0,-5,'S','FontSize',14,'HorizontalAlignment','Center');
	HT2=text(90,-5,'W','FontSize',14,'HorizontalAlignment','Center');
   HT3=text(180,-5,'N','FontSize',14,'HorizontalAlignment','Center');
   HT4=0;
elseif (isequal(arg,'N') == 1)
   % North
   set(gca,'Position',[0.05 0.07 0.90 0.55],'Color',backCol);
   axis([90 270 0 70]);
	HT1=text(90,-5,'W','FontSize',14,'HorizontalAlignment','Center');
	HT2=text(180,-5,'N','FontSize',14,'HorizontalAlignment','Center');
   HT3=text(270,-5,'E','FontSize',14,'HorizontalAlignment','Center');
   HT4=0;
elseif (isequal(arg,'O') == 1)
   % East
   set(gca,'Position',[0.05 0.07 0.90 0.55],'Color',backCol);
   axis([180 360 0 70]);
	HT1=text(180,-5,'N','FontSize',14,'HorizontalAlignment','Center');
	HT2=text(270,-5,'E','FontSize',14,'HorizontalAlignment','Center');
	HT3=text(360,-5,'S','FontSize',14,'HorizontalAlignment','Center');
   HT4=0;
elseif (isequal(arg,'S') == 1)
   % South
   set(gca,'Position',[0.05 0.07 0.90 0.55],'Color',backCol);
   axis([-90 90 0 70]);
	HT1=text(-90,-5,'E','FontSize',14,'HorizontalAlignment','Center');
	HT2=text(0,-5,'S','FontSize',14,'HorizontalAlignment','Center');
	HT3=text(90,-5,'W','FontSize',14,'HorizontalAlignment','Center');
   HT4=0;
   negindC = find(azC >= 270);
   azC(negindC) = azC(negindC)-360;
   negind = find(az >= 270);
   az(negind) = az(negind)-360;
   negindP = find(azP >= 270);
   azP(negindP) = azP(negindP)-360;
   if (azM >= 270)
      azM = azM-360;
   end
end

hold on

% plot names of constellations
if (showconst == 1)
   for c = 1:anz_bilder
      text(azC(c),hC(c),sternbilder(c).name,'Color','c','HorizontalAlignment','Center','FontSize',7);
   end
end

% plot stars by magnitude
plot(az(ind00),h(ind00),'w.','MarkerSize',19)
plot(az(ind01),h(ind01),'w.','MarkerSize',14)
plot(az(ind02),h(ind02),'w.','MarkerSize',9)
plot(az(ind03),h(ind03),'w.','MarkerSize',4)

% text alignment
HorAlign = horalign([azP azM],[hP hM],10.0);

% plot planets incl. sun
for p = 9:-1:1
   if (p == 3)
      plot(azP(p),hP(p),'y.','MarkerSize',24);
   else
      plot(azP(p),hP(p),'m.','MarkerSize',19);
	end
   if (HorAlign(p) == 1)
      HA = 'left';
   else
      HA = 'right';
   end
   text(azP(p),hP(p),planetenStr(p,:),'Color',[1 1 1],'HorizontalAlignment',HA);
end

% plot moon
plot(azM,hM,'y.','MarkerSize',24);
if (HorAlign(10) == 1)
   HA = 'left';
else
   HA = 'right';
end
text(azM,hM,' MOON ','Color',[1 1 1],'HorizontalAlignment',HA);

% Enable Pick
set(H_PICK,'Enable','on');

Contact us at files@mathworks.com