Чего получил-то? Сам-то хоть понял?
да вот:
%-------------------------------------------------------------------------------------------------------------
% Drawing a horoscope of a date.
% Input : the date of a horoscope
% Reference:
http://naif.jpl.nasa.gov/pub/naif/toolkit_docs/MATLAB/% Reference:
http://ru.scribd.com/doc/57253920/Anatoly-T-fomenko-Mysteries-of-Egyptian-Zodiacs-and-Other-Riddles-of-Ancient-History% Tested : Matlab R2011b, 24.04.2013 - 09.07.2013
% By : Poltavsky Sasha
% URL :
http://www.bible-exodus.narod2.ru/articles/astro_ephemeris/jpl_ephemeris/jpl_ephemeris.html#zodiacs%--------------------------------------------------------------------------------------------------------------
% ATTENTION !!! ВНИМАНИЕ !!!
% before runing tha script download and install additional files : % перед тем как запустить скрипт, скачайте и добавьте
%
ftp://ssd.jpl.nasa.gov/pub/eph/planets/bsp/de422.bsp - ephemeris for 2999 BC ... 3000 AD
%
ftp://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/a_old_versions/de408.bsp - ephemeris for 10000 BC ... 10000 AD
%
http://www.bible-exodus.narod2.ru/articles/astro_ephemeris/jpl_ephemeris/data/my_sites.tpc - or use '399' instead of 'observer'
%
http://www.bible-exodus.narod2.ru/articles/astro_ephemeris/jpl_ephemeris/data/my_frames.tf - or use '399' instead of 'observer'
%
http://naif.jpl.nasa.gov/pub/naif/pds/data/ros-e_m_a_c-spice-6-v1.0/rossp_1000/DATA/FK/RSSD0002.TF%
http://bible-exodus.narod2.ru/articles/astro_ephemeris/jpl_ephemeris/scripts/mice_toolkit_n0064/mice_tol_1e-4_lib.zip - for increasing TOL
% add paths to standart.tm
% also you need to add my subroutines : set_observer.m, get_observer.m,
% angl2minsec.m ... - they are below in commentc (uncomment Ctrl + T, cut,
% and save with it's file names, add their paths to MatLab
clear all;
cspice_furnsh( 'd:\Documents and Settings\user\Мои документы\MATLAB\astro\mice\data\standard.tm' );
fl_lng = 1; % Output language English = 1, Russian =0, Both = 2
% the horoscope epoch % время гороскопа
et_horos = cspice_str2et ( '18:07 12-23-0001 BC' );
% et_horos = cspice_str2et ( datestr(now, 'mmmm dd, yyyy HH:MM:SS UTC+2'));
% et_horos = cspice_str2et ( '2013-07-03 23:59 UTC+0'); % horoscope epoch - 12-23-0001 BC = Xmas % момент гороскопа
% the epoch for scale of second zodiac (for publishing date, for example
et_zodiac_2=et_horos; % ввести свое значние et_zodiac_2 для его отрисовки
%et_zodiac_2 = cspice_str2et ( '1185-03-27 06:00 UTC+2'); % publication epoch % момент пересчета гороскопа
% set observer coordinates or , for example 'CAIRO' settled in my_sites.bsp % место наблюдения
[observer, topo_frame ]=set_observer(44.6,3.53333,0); % Севастополь (близ мыса Фиолент) - для гороскопа Рождества Христова по ФиН
% [observer, topo_frame ]=set_observer(45,0,0); % широта, долгота, высота - LAN, LON, ALT
[lat, lon, alt, sslat, sslon] = get_observer (observer); % for local solar time % узнать долготу места для местного солн.времени
% to set the time format % установить формат времени (JCAL - Julian calendar % Юлианский календарь)
TIMFMT = 'DD.MM.YYYY ERA HR:MN ::RND::JCAL::UTC'; TIMLEN = 35;
abcr = 'LT+S'; step = cspice_spd/2; crds='CYLINDRICAL' ; crd = 'LONGITUDE' ; adj = 0.0; MAXWIN = 20000; MAXIVL = MAXWIN / 2;
% % adjustment on precession = 360° for 25776 years % поправки на прецессию 360° за 25776 лет
% % used 1976 IAU precession model, built into SPICE % использована модель "1976 IAU" встроенная в SPICE.
% [pos1, ~] = cspice_spkpos ( 'SUN', et_horos, 'ECLIPDATE' , abcr, observer); % позиция планеты от эклиптики той эпохи
% [pos2, ~] = cspice_spkpos ( 'SUN', et_horos, 'ECLIPJ2000', abcr, observer); % позиция планеты по нынешней эклиптике
% [~, lon1, ~] = cspice_recrad(pos1); [~, lon2, ~] = cspice_recrad(pos2); lon3 = lon2-lon1;
% prcs_at_horos = cspice_convrt(lon3,'radians','degrees');
% [pos1, ~] = cspice_spkpos ( 'SUN', et_zodiac_2, 'ECLIPDATE' , abcr, observer); % позиция планеты от эклиптики той эпохи
% [pos2, ~] = cspice_spkpos ( 'SUN', et_zodiac_2, 'ECLIPJ2000', abcr, observer); % позиция планеты по нынешней эклиптике
% [~, lon1, ~] = cspice_recrad(pos1); [~, lon2, ~] = cspice_recrad(pos2); lon3 = lon2-lon1;
% prcs_zodiac_2 = cspice_convrt(lon3,'radians','degrees');
% the same easy way % более просто способ
prcs_at_horos=360/((cspice_tyear*25776)/(now-et_horos)); prcs_zodiac_2=360/((cspice_tyear*25776)/(now-et_zodiac_2));
% data array for drawing the zodiac % данные для отрисовки зодиака
planets = {'Jup', 'Sat', 'Ven', 'Mars', 'Mer', 'Sun', 'Moon', 'Earth'; ...
'JUPITER BARYCENTER', 'SATURN BARYCENTER', 'VENUS', 'MARS', 'MERCURY', 'SUN', 'MOON', 'EARTH'};
dots = [15,12,10,10,10,35,35,1];
colors = [[0.2 0 0]; [0 0.2 0.3]; [0 0.2 0.1]; [0.6 0 0]; [0.4 0 0]; [0.7 0.5 0]; [0.5 0.4 0.1]; [0.8 0.4 0.1]];
lats = [0,0,0,0,0,0,0,0]; lats_h = [0,0,0,0,0,0,0,0]; lons = [0,0,0,0,0,0,0,0];
sky_s = {'темно - ночь', 'астрономические сумерки', 'навигационные сумерки', 'гражданские сумерки', 'светло - день'};
sky_e = {'dark - night', 'astronomical twilight', 'nautical twilight', 'civil twilight', 'daylight'};
zodiac_names = {'<-Овен ', 'Телец ', 'Близнецы ', 'Рак ', 'Лев ', ...
'Дева ', 'Весы ', 'Скорпион ', 'Стрелец ', 'Козерог ', 'Водолей ', 'Рыбы ';
'<-Aries', 'Taurus', 'Gemini', 'Cancer', 'Leo', 'Virgo', 'Libra', ...
'Scorpio', 'Sagittarius', 'Capricorn', 'Aquarius', 'Pisces';
'<-Ari|Овн', 'Tau|Тлц', 'Gem|Блз', 'Cnc|Рак', 'Leo|Лев', 'Vir|Дева', ...
'Lib|Весы', 'Sco|Скрпн', 'Sgr|Стрлц', 'Cap|Кзрг', 'Aqr|Вдл', 'Psc|Рыбы'};
% Ecliptics coordinates of the East point of the local topo frame - восток в координатах эклиптики
mat1 = cspice_pxform( topo_frame, 'ECLIPJ2000', et_horos );
mat2 = mat1 * cspice_rotate( -0.5*cspice_pi, 3); mat2 = mat2(:,1,:);
[~, lon_east, lat_east] = cspice_reclat(mat2);
lon_east = 180 + lon_east * cspice_dpr; lat_east = - lat_east * cspice_dpr;
% Ecliptics coordinates of the West point of the local topo frame - запад в координатах эклиптики
mat2 = mat1 * cspice_rotate( 0.5*cspice_pi, 3); mat2 = mat2(:,1,:);
[~, lon_west, lat_west] = cspice_reclat(mat2);
lon_west = 180 + lon_west * cspice_dpr; lat_west = - lat_west * cspice_dpr;
% Ecliptics to horizon inclination - наклон эклиптики к горизонту
mat1 = cspice_pxform( topo_frame, 'ECLIPJ2000', et_horos ); mat1 = mat1(:,3,:);
[~, lon_p, lat_p] = cspice_reclat(mat1);
eclip_horiz_incl = - 90 + lat_p * cspice_dpr;
% North pole longitude in ecliptic
north_pole_lon_e = 180 + lon_p * cspice_dpr;
% Longitudes of intercection of ecliptic and horizon
lon_w_e = north_pole_lon_e - 90 ; if lon_w_e < 0 ; lon_w_e = lon_w_e + 360; end
lon_e_e = north_pole_lon_e + 90 ; if lon_e_e > 360 ; lon_e_e = lon_e_e - 360; end
% разница в азимуте между точками пересечения эклиптики и горизонта к точкам востока и запада
d_lon = asind( tand (10) / tand (eclip_horiz_incl) );
% output an image of goroscope % вывод гороскопа
p_place = horzcat(observer, sslat, sslon, ' (', angl2minsec(lat,2), ', ', angl2minsec(lon,2), ') ');
p_epoch = horzcat( cspice_timout( et_horos, 'DD.MM.YYYY ERA JCAL HR:MN ::RND::JCAL::UTC UTC '),...
' ( ', cspice_timout(et_horos, 'DD.MM.YYYY ::UTC'),' ) ', cspice_et2utc( et_horos, 'J', 2 ));
if fl_lng == 0
p_name_t=horzcat('Планеты в созвездиях на ', p_epoch, ' место ', p_place);
elseif fl_lng > 0
p_name_t=horzcat('Planets in constellations at ',p_epoch,' in ', p_place);
end
p_name={p_name_t;' '};
scrsz = get(0,'ScreenSize');
h = figure('Name', p_name_t ,'NumberTitle','off', 'Position',[1 scrsz(4)/2.6 scrsz(3) (scrsz(4)/2)-20]);
% цвет неба в зависимсти от наличия солнца на небе, the color of sky relatively sun's presens
[pos, ~] = cspice_spkpos ( 'SUN', et_horos, topo_frame, abcr, observer);
[~, ~, sun_alt] = cspice_reclat(pos); sun_alt=sun_alt * cspice_dpr;
if fl_lng == 0; sky_t = horzcat('Высота солнца',angl2minsec(sun_alt,0),', ');
else sky_t = horzcat('Altitude of sun is',angl2minsec(sun_alt,0),', '); sky_s = sky_e;
end
if sun_alt < -12 && sun_alt > -18; sky_color= [0.5 0.4 0.6]; sky_s = horzcat(sky_t, cell2mat(sky_s(2))) ;
elseif sun_alt < -7 && sun_alt > -12 ; sky_color= [0.8 0.7 0.7]; sky_s = horzcat(sky_t, cell2mat(sky_s(3))) ; % навигационные сумерки
elseif sun_alt < 0 && sun_alt > -7 ; sky_color= [0.9 0.9 0.9]; sky_s = horzcat(sky_t, cell2mat(sky_s(4))) ; % гражданские сумерки
elseif sun_alt > 0; sky_color= [1 1 1]; sky_s = horzcat(sky_t, cell2mat(sky_s(5))) ; % светло - день
elseif sun_alt < -18 ; sky_color= [0.3 0.4 0.6]; sky_s = horzcat(sky_t, cell2mat(sky_s(1))) ;
end % темно - ночь
% % drawing the Earth's body - тело цельное, слева восток - не касается границы графика, справа - запад
% if (lon_e_e + d_lon) < 360 && (lon_e_e - d_lon) < 360 && (lon_w_e + d_lon) > 0 && (lon_w_e - d_lon) > 0 && (lon_e_e > lon_w_e)
e1=lon_e_e + d_lon ; e2 = lon_e_e - d_lon ; w1=lon_w_e + d_lon ; w2 = lon_w_e - d_lon ; y=[10 -10 -10 10 ]; y2=y; y3=y;
if lon_e_e > lon_w_e
x=[e1 e2 w1 w2]; x2=x()+360; x3=x()-360;
elseif lon_e_e < lon_w_e
x=[e1 e2 w1-360 w2-360]; x2= [e1+360 e2+360 w1 w2]; x3=[]; y3=[];
end
fill([360 360 0 0 ], y, sky_color); hold on
f1=fill(x, y, [.7 .8 .7], x2, y2, [.7 .8 .7], x3, y3, [.7 .8 .7]); set(f1,'facealpha',.5)
title(p_name);
% calculating of positions of planets % расчет позиций планет
if fl_lng == 0
fprintf ( '\n\n позиции планет в эклиптике J2000.0, на %s \n\n', cspice_timout( et_horos, 'DD.MM.YYYY ERA JCAL HR:MN ::RND::JCAL::UTC UTC')) ;
else
fprintf ( '\n\n the positions of planet at ecliptics J2000.0, на %s \n\n', cspice_timout( et_horos, 'DD.MM.YYYY ERA JCAL HR:MN ::RND::JCAL::UTC UTC')) ;
end
fl=0; %
n=numel(planets(1,:)); % количество записей
for ii=1:n-1
target=planets(2,ii);
[pos, ~] = cspice_spkpos ( target, et_horos, 'ECLIPJ2000', abcr, observer); % позиция планеты
[~, lon_p, lat_p] = cspice_recrad(pos); % картезианские координаты в радиальные небесной сферы
lons(ii) = lon_p * cspice_dpr; lats(ii)= lat_p * cspice_dpr ;
[pos, ~] = cspice_spkpos ( target, et_horos, topo_frame, abcr, observer);
[~, ~, lat_h] = cspice_reclat(pos); lats_h(ii) =lat_h * cspice_dpr; %az=-(-180+lon);
fprintf ( '%s LON %s\t LAT %s\t \n', cell2mat(target), num2str(lons(ii)), num2str(lats(ii)));
text(360,-18- ii*5, planets(1,ii)); text(345,-18- ii*5, horzcat('Lon ', num2str(lons(ii))));
text(300,-18- ii*5, horzcat('Lat ', num2str(lats(ii))));
text(260,-18- ii*5, horzcat('Alt ', num2str(lats_h(ii))));
end
% drawing planets % отрисовка планет
n=numel(lons(:)); % количество записей
for ii=1:n-1
text(lons(ii)+1,lats(ii)+3, planets(1,ii), 'HorizontalAlignment','right')
p=scatter(lons(ii), lats(ii), dots(ii), colors(ii,:), 'filled');
end
if fl_lng == 0
text(360,-65,'Шкалы: красным начала созвездий, синим границы зодиака на момент гороскопа. Спасибо НАСА! ® Полтавский Саша, vasnas@live.com');
text(180, -25, sky_s);
else
text(360,-65,'Scales: red - begining of constellations, blue - zodiac''s segments at epoch of horoscope. Thanks to NASA! ® Poltavsky Sasha, vasnas@live.com');
text(180, -25, sky_s);
end
ax1=gca; axis image ; set(ax1,'XDir','reverse', 'Color', [0.9 0.9 1],'XColor', [0.5 0 0], 'Layer', 'top');
% beginings of constellations in ecliptics at J2000.0 % начала созвездий по эклиптике J2000
set(ax1,'Xtick', [31, 56, 92, 118, 137, 172, 215, 236, 266, 296, 326, 349], 'XLim',[0,360], 'YLim',[-10,10]);
% names of constellations % названия создвездий
if fl_lng==0 ; zodiac_names_out = zodiac_names(1,:);
elseif fl_lng == 1 ; zodiac_names_out = zodiac_names(2,:);
elseif fl_lng == 2; zodiac_names_out = zodiac_names(3,:);
end
set(gca,'XTickLabel',zodiac_names_out)
grid on
% scale of zodiac at goroscope epoch - blue % шкала зодиака на эпоху гороскопа - синяя
Ax2=axes('Position', get(ax1,'Position'),'XAxisLocation','top'); axis image ;
set(Ax2,'XLim',[0 360], 'YLim',[-10 10], 'color','none','XColor', [0 0 0.7]) % 'xtick',(prcs_at_horos:30:360+prcs_at_horos)
set(Ax2, 'xticklabel', zodiac_names_out,'xtick', prcs_at_horos:30:360+prcs_at_horos,'XDir','reverse','FontSize',8 ,'Yticklabel','')
grid on
if et_horos ~= et_zodiac_2
% scale of zodiac at publications (or anything else) epoch - green % шкала зодиака на эпоху публикации - зеленая
Ax3=axes('Position',get(ax1,'Position'),'XAxisLocation','top'); axis('image'); set(Ax3,'XLim',[0 360], 'YLim',[-10 10]);
set(Ax3,'color','none'); set(Ax3,'XColor', [0 0.7 0])
set(Ax3,'xtick',(prcs_zodiac_2:30:360+prcs_zodiac_2),'xticklabel',(0:30:360),'XDir','reverse');
set(Ax3,'Yticklabel','','TickDir', 'out','FontSize',8)
grid on
end
% scale of zodiac at publications (or anything else) epoch - green % шкала зодиака на эпоху публикации - зеленая
Ax4=axes('Position',get(ax1,'Position'),'XAxisLocation','bottom'); axis('image'); set(Ax4,'XLim',[0 360], 'YLim',[-10 10]);
set(Ax4,'color','none','XColor', [0 0.2 0])
set(Ax4,'xtick', (0:30:360),'xticklabel',(0:30:360),'XDir','reverse');
set(Ax4,'Yticklabel','','TickDir', 'out','FontSize',7)
grid on
box off
cspice_kclear
clear all
% %-------------------------------------------------------------------------------------------------------------
% % Getting coordinates of observer, after function set_observer(). Output in degree,
% % Reference:
http://naif.jpl.nasa.gov/pub/naif/toolkit_docs/MATLAB/% % Tested : Matlab R2011b
% % By : Poltavsky Sasha 04.05.2012
% % URL :
www.bible-exodus.narod2.ru/articles/astro_ephemeris/jpl_ephemeris/jpl_ephemeris.html#get_kernel_data% %--------------------------------------------------------------------------------------------------------------
%
% function [lat, lon, alt, sslat, sslon] = get_observer (observer, varargin)
% if nargin == 0 ; observer = 'DYN'; end
%
% [val, found] = cspice_gdpool( horzcat('TKFRAME_',observer, '_TOPO_ANGLES'), 2, 1 ); lat=90+val;
% if ~found; error(horzcat('the latitude of ',observer,' was not found in the pool!')); end
% [val, found] = cspice_gdpool( horzcat('TKFRAME_',observer, '_TOPO_ANGLES'), 1 ,1); lon= - val;
% if ~found; error(horzcat('the longtitude of ',observer,' was not found in the pool!')); end
% %[val, found] = cspice_gdpool( horzcat('TKFRAME_',observer, '_TOPO_ALT'), 1 ,1); alt = val;
% %[alt, found] = cspice_gdpool( horzcat(observer, '_LATLON'), 3 ,1);
% %if ~found; error(horzcat('the altitude of ',observer,' was not found in the pool!')); end
% alt=0;
%
%
%
% if lat < 0; sslat = sprintf('% 2.2fS',abs(lat)); else sslat = sprintf('% 2.2fN',abs(lat)); end
% if lon < 0; sslon = sprintf('% 2.2fW',abs(lon)); else sslon = sprintf('% 2.2fE',abs(lon)); end
% end
%
% %-------------------------------------------------------------------------------------------------------------
% % Setting coordinates of observer, input in degree,
% % Reference:
http://naif.jpl.nasa.gov/pub/naif/toolkit_docs/MATLAB/% % Tested : Matlab R2011b
% % By : Poltavsky Sasha 04.05.2012
% % URL :
www.bible-exodus.narod2.ru/articles/astro_ephemeris/jpl_ephemeris/jpl_ephemeris.html#get_kernel_data% %--------------------------------------------------------------------------------------------------------------
% function [observer, frame] = set_observer (lat, lon, alt, varargin)
% if nargin == 1; name = lat;
% [~ , found] = cspice_bods2c(name);
% if (found); observer = name; frame = strcat(observer,'_TOPO');
% else
% disp('The point of surface or the frame DYN_TOPO for the body are not specified')
% end
% else
% if nargin == 2; observer = 'DYN'; frame = strcat(observer,'_TOPO'); alt=0; end
% if nargin == 3; observer = 'DYN'; frame = strcat(observer,'_TOPO'); end
%
% cspice_pdpool('TKFRAME_DYN_TOPO_ANGLES',[-lon; -90 + lat; 180]);
% cspice_lmpool(horzcat('DYN_LATLON = (', num2str(lat), ',', num2str(lon), ',', num2str(alt), ' )' ));
% %cspice_lmpool(horzcat('TKFRAME_DYN_TOPO_ALT = (',num2str(alt),' )'));
%
%
%
% end
%
%
% return
% %-------------------------------------------------------------------------------------------------------------
% % Function of translating decimal degrees into sexagesimal system
% % Tested : Matlab R2011b
% % By : Poltavsky Sasha 02.06.2012
% % URL :
www.bible-exodus.narod2.ru/articles/astro_ephemeris/jpl_ephemeris/jpl_ephemeris.html#angle_10_to_60% %--------------------------------------------------------------------------------------------------------------
% % an - angle in decimal degrees
% % prec - precision in "
%
% function anstr = angl2minsec (an, prec, varargin)
%
% if nargin == 1 ; sprec='0';
% else
% sprec=int2str(prec);
% end
%
% if an<0 ; s_sign='-'; else s_sign=' '; end
%
% an1=fix(abs(an));
% an2=abs(rem(an,1));
% an2=an2*0.6;
% an3=rem(an2*100,1)*100;
% an2=fix(an2*100);
% an3=an3*0.6;
% anstr=sprintf('%s%s°%s''%s"', s_sign, int2str(an1), int2str(an2), num2str(an3, strcat('%2.',sprec,'f')));
% return
