%------------------------------------------------------------------------------------------------------------- % 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.4 0 0]; [0 0.4 0.6]; [0.1 0.5 0.3]; [0.6 0 0]; [0.8 0 0]; [0.7 0.5 0]; [0.8 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