<< 10. Приложение D | Оглавление | Литература >>
Разделы
11. Приложение E
11.1. Тексты программ на языке FORTRAN 90
Листинг E 3.1. Модуль чтения каталога Hipparcos
MODULE HipMain IMPLICIT NONE ! Расположение полной версии каталога Hipparcos CHARACTER(*), PARAMETER :: HipparcosName ='D:/HIP/hip\_main.dat' INTEGER, PARAMETER :: HipNumOfStars = 118218 ! Число звезд INTEGER, PARAMETER :: u = 10 ! Номер файла TYPE THipparcos SEQUENCE INTEGER(4) :: HIP ! Номер звезды по Hipparocs ! Астрометрическая информация REAL(8) :: RAdeg,DEdeg ! экваториальные координаты в градусах REAL(8) :: Plx ! тригонометрический параллакс в mas REAL(8) :: pmRa,pmDE ! собственные движения ma*cos(d) и md CHARACTER(1) :: AstroRef ! Флаг для кратных систем ! Фотометрическая информация REAL(4) :: VMag ! Звездная вел. по шкале Джонсона REAL(4) :: B\_V ! Показатель цвета B-V по шкале Джонсона ! Ошибки соответствующих величин REAL(8) :: sigma\_RAdeg,sigma\_DEdeg REAL(8) :: sigma\_Plx REAL(8) :: sigma\_pmRa,sigma\_pmDE CHARACTER(10) Sp ! Развернутый спектральный класс LOGICAL NoRaDe ! Нет данных о точных координатах LOGICAL NoPlx ! Нет данных о параллаксе LOGICAL Nopm ! Нет данных о собственных движениях LOGICAL NoVMag ! Нет данных о звездной величине LOGICAL NoB\_V ! Нет данных о показателе цвета END TYPE THipparcos CONTAINS SUBROUTINE OpenHipparcosMain ! Открытие файла каталога OPEN(u, file = HipparcosName) END SUBROUTINE OpenHipparcosMain SUBROUTINE CloseHipparcosMain ! Закрытие файла каталога CLOSE(u) END SUBROUTINE CloseHipparcosMain LOGICAL FUNCTION ReadHipparcosMain(s) ! Чтение данных о звезде TYPE(THipparcos), INTENT(out) :: s CHARACTER(450) hs ! Запись строки каталога IF (EOF(u)) THEN ReadHipparcosMain=.false. RETURN ELSE ReadHipparcosMain=.true. END IF READ(u,'(A450)') hs ! Чтение одной строки каталога ! Сбрасываем флаги событий s.NoRaDe = .False. s.NoPlx =.False. s.Nopm =.False. s.NoVMag =.False. s.NoB\_V =.False. ! Интерпретация с 12 байт, начиная с 3-го - это номер HIP read(hs(3:14),*) s.hip ! Чтение координат: по 12 байт с 52 и с 65 позиции ! Функция TRIM удаляет из строки пробелы, а LEN возвращает длину ! строки, соответственно, если это 0, то в строке только пробелы IF (LEN(TRIM(hs(52:63))) == 0) THEN s.NoRaDe=.true. s.RADeg=0.0 ! на всякий случай записываем 0 ELSE READ(hs(52:63),*) s.RAdeg END IF IF (LEN(TRIM(hs(65:76))) == 0) THEN s.NoRaDe=.true. s.DEDeg=0.0 ELSE read(hs(65:76),*) s.DEdeg END IF ! Чтение параллакса - 7 байт с 80-й позиции IF (LEN(TRIM(hs(80:86))) == 0) THEN s.NoPlx=.true. s.Plx=0.0 ELSE read(hs(80:86),*) s.Plx END IF ! Чтение собственных движений: по 8 байт с 88 и с 97 позиции IF (LEN(TRIM(hs(88:95))) == 0) THEN s.NoPM=.true. s.pmRA=0.0 ELSE read(hs(80:86),*) s.pmRA END IF IF (LEN(TRIM(hs(88:95))) == 0) THEN s.NoPM=.true. s.pmDE=0.0 ELSE read(hs(97:104),*) s.pmDE END IF s.AstroRef=hs(78:78) ! Флаг кратной звезды ! Чтение зв.величины и показателя цвета B-V по шкале Джонсона IF (LEN(TRIM(hs(42:46))) == 0) THEN s.NoVMag=.true. s.VMag=0.0 ELSE read(hs(42:46),*) s.VMag END IF IF (LEN(TRIM(hs(246:251))) == 0) THEN s.NoB\_V=.true. s.B\_V=0.0 ELSE read(hs(246:251),*) s.B\_V END IF ! Данные об ошибках всегда есть, если присутствуют сами величины IF (.NOT. s.NoRADE) THEN read(hs(106:111),*) s.sigma\_RAdeg read(hs(113:118),*) s.sigma\_DEdeg ENDIF IF (.NOT. s.NoPlx) THEN read(hs(120:125),*) s.sigma\_Plx END IF IF (.NOT. s.Nopm) THEN read(hs(127:132),*) s.sigma\_pmRA read(hs(134:139),*) s.sigma\_pmDE END IF s.Sp=hs(436:445) ! Чтение данных о спектральном классе END FUNCTION ReadHipparcosMain END MODULE HipMain
Листинг E 3.2. Подсчет звезд без данных о координатах, собственных движениях и параллаксах
PROGRAM Test2 USE HipMain INTEGER(4) :: NoCoord = 0 ! Счетчик звезд без точных координат INTEGER(4) :: NoProp = 0 ! Счетчик звезд без собств. движений INTEGER(4) :: NoPar = 0 ! Счетчик звезд без параллаксов TYPE(THipparcos) :: s CALL OpenHipparcosMain DO WHILE (ReadHipparcosMain(s)) ! Сравнение логических переменных IF (s.NoRADE) NoCoord= NoCoord+1 IF (s.Nopm) NoProp = NoProp+1 IF (s.NoPlx) NoPar = NoPar+1 END DO CALL CloseHipparcosMain print *,'No coord ',NoCoord print *,'No PM ',NoProp print *,'No Par ',NoPar read *,i END PROGRAM Test2
Листинг E 3.3. Вычисление распределения звезд по абсолютной звездной величине
PROGRAM AbsMagDistrib USE HipMain INTEGER, PARAMETER :: LOW = -12, HIGH=+12 TYPE(THipparcos) :: s; INTEGER(4) :: a(-12:+12) ! статистика INTEGER(4) :: i ! вспомогательная переменная REAL(8) :: r ! расстояние REAL(8) :: m ! абсолютная звездная величина FORALL (i=-12:+12) A(i)=0 CALL OpenHipparcosMain(); DO WHILE (ReadHipparcosMain(s)) IF (s.NoPlx) CYCLE ! нет данных о паралл. IF (s.Plx<=0.0) CYCLE ! неположительный параллакс IF (s.sigma\_plx/s.plx>0.5) CYCLE ! точность хуже 50\% r=1000.0/s.plx; ! Вычисление расстояния в пк m=S.VMag-5.0*log10(r)+5.0 ! Вычисл. абсолютной звезд. величины i=FLOOR(m+0.5) ! Определение индекса ячейки массива IF ( (i>=low) .and. (i<=high)) a(i)=a(i)+1 ! ув.на 1 END DO CALL CloseHipparcosMain() PRINT '(I3,1X,I7)',(i,a(i),i=-12,12) READ *,i END PROGRAM
Листинг E 3.4. Вычисление средней абсолютной звездной величины звезд, список которых находится в файле lumin.txt
PROGRAM AVER Use HipMain TYPE(THipparcos) s REAL(8) :: r ! расстояние REAL(8) :: m ! абсолютная звездная величина REAL(8) :: mav = 0.0 ! средняя абсолютная звездная величина INTEGER(4) :: n = 0 ! количество подходящих звезд CALL OpenHipparcosMain ! Инициализация критерия print '(I6 " звезд в критерии.")',InitCriteria('lumin.txt') DO WHILE (ReadHipparcosMain(s)) IF (s.NoPlx) CYCLE ! нет данных о паралл. IF (s.plx<=0.0) CYCLE ! неположительный параллакс r=1000.0/s.plx ! Вычисление расстояния в пк m=S.VMag-5.0*log10(r)+5.0 ! Вычисл. абс. звезд. величины IF (inCelestia(s.HIP)) THEN ! Звезда в списке Celestia mav=mav+m ! накопление суммы абс. зв. величин n=n+1 ! суммирование числа звезд END IF END DO CALL ClearCriteria ! Очистка критерия CALL CloseHipparcosMain mav=mav/n ! Вычисление среднего значения print '("Средняя абсолютная звездная величина ",F6.2,".")',mav print '("Обработано ",I6," звезд.")',n END PROGRAM
Листинг E 3.5. Исходный текст функций InitCriteria, InCelestia, ClearCriteria
Добавления в раздел описаний модуля HipMain
INTEGER(4), ALLOCATABLE :: List(:) INTEGER(4) :: NList = 0
Добавления в раздел процедур модуля HipMain
INTEGER(4) FUNCTION InitCriteria(name) CHARACTER (*) :: name INTEGER, PARAMETER :: t = 11 ! файл INTEGER(4) :: i ! индекс для цикла do OPEN(t, file = name, mode='read') ! Открываем файл DO i=1,2 ! Пропуск двух первых строк READ(t,*) END DO READ(t,*) NList ! В третьей строке - количество звезд PRINT *,NList ALLOCATE(List(NList)) ! Выделение памяти под список DO i=4,12 ! Пропуск с 4 по 12 строку READ(t,*) END DO DO i=1,NList ! Чтение номеров звезд READ(t,'(2X,I12)') List(i) END DO CLOSE(t) InitCriteria=NList END FUNCTION InitCriteria ! Очистка критерия SUBROUTINE ClearCriteria DEALLOCATE(List) NList=0 END SUBROUTINE ClearCriteria ! Функция проверяет, есть ли звезда в списке LOGICAL FUNCTION InCelestia(n) INTEGER(4) :: n INTEGER(4) :: i InCelestia=.false. ! объект пока не найден IF (NList==0) return ! если критерий не установлен - выход DO i=1,NList ! обход звезд в цикле do IF (List(i)==n) THEN ! если звезда найдена в списке InCelestia=.true. EXIT ! досрочно прервать цикл END IF END DO END FUNCTION InCelestia
Листинги E 4.1-4.3. Процедуры перевода координат и отрисовки координатной сетки (Перевод в экранные координаты не требуется, поскольку в стандартной библиотеке фортрана можно выбрать любую удобную декартову систему с вещественными значениями координат)
MODULE Projection USE DFLIB IMPLICIT NONE REAL(8), PARAMETER :: PI = 3.1415926535897932384626433832795 CONTAINS SUBROUTINE Aitoff(l,b,x,y) REAL(8), INTENT(IN) :: l,b ! Сферические координаты в радианах REAL(8), INTENT(OUT) :: x,y ! Декартовы координаты REAL(8) :: s, l1 ! Вспомогательные переменные IF (l>PI) THEN ! Приведение l в диапазон -Pi до +Pi l1=l-2*Pi ELSE l1=l END IF S=sqrt(1.0+cos(b)*cos(l1/2)) ! Знаменатель формул (4.1) x=-2*cos(b)*sin(l1/2)/s y=sin(b)/s END SUBROUTINE Aitoff REAL(8) FUNCTION radi(x) ! Перевод градусов в радианы INTEGER, INTENT(IN) :: x radi=x/180.0*PI END FUNCTION radi REAL(8) FUNCTION rad(x) ! Перевод градусов в радианы REAL(8), INTENT(IN) :: x rad=x/180.0*PI END FUNCTION rad SUBROUTINE AitoffGrid (Step,Gr) INTEGER, INTENT(IN) :: Step ! Шаг сетки в градусах LOGICAL, INTENT(IN) :: Gr ! Флаг - в градусах или в часах INTEGER :: i,j ! Переменные циклов do REAL(8) :: l,b ! Галактические координаты REAL(8) :: x,y ! Декартовы координаты CHARACTER(8) s ! Строка для подписей INTEGER :: h ! Для разметки осей TYPE (wxycoord) wxy INTEGER(2) :: status2 INTEGER(4) :: status4 ! Нанесение сетки меридианов status4 = SetColorRGB({\#00FF00) DO i=-180,+180,Step l=radi(i) ! Перевод в радианы DO j=-90,+90,5 ! Цикл построения вдоль меридиана ! Вычисление точки меридиана b=radi(j) ! Перевод в радианы широты CALL Aitoff(l,b,x,y) ! Перевод в декартовы координаты ! Если точка первая (j=-90), то помещаем графический курсор ! в точку (x,y) функцией MoveTo\_W, если точка не первая, то ! "прочерчиваем" курсором линию из предыдущей точки ! в точку (u,v) функцией LineTo\_W. IF (j==-90) THEN CALL MoveTo\_W(x,y,wxy) ELSE status2=LineTo\_W(x,y); END IF END DO ! j END DO ! i ! Нанесение сетки параллелей - аналогично предыдущему DO j=-90,+90,Step b=radi(j) DO i=-180,+180,5 ! цикл построения вдоль параллели l=radi(i) CALL Aitoff(l,b,x,y); IF (i==-180) THEN CALL MoveTo\_W(x,y,wxy) ELSE status2=LineTo\_W(x,y) END IF END DO ! i END DO ! j status2=SetFont('t"Arial"h10') status4 = SetColorRGB({\#FFFFFF) ! Подписи меридианов вдоль экватора DO i=-180,+180,Step ! Вычисление координаты точки вывода надписи l=Radi(i); CALL Aitoff(l,0.0\_8,x,y) ! Если Gr истина, то разметка в градусах, иначе - в часах IF (Gr) THEN h=i ELSE h=i/15 IF (h<0) h=h+24; END IF write(s,'(I4)') h ! Преобразование значения h в текст. строку Call MoveTo\_W(x, y, wxy) Call OUTGTEXT(s) END DO ! Подписи параллелей вдоль нулевого меридиана - аналогично DO j=-90,+90,Step IF (j /= 0) THEN ! Экватор не подписываем b=Radi(j); CALL Aitoff(0.0\_8,b,x,y) write(s,'(I4)') j Call MoveTo\_W(x, y, wxy) Call OUTGTEXT(s) END IF END DO END SUBROUTINE AitoffGrid SUBROUTINE Galaxy(a,d,l,b) REAL(8), INTENT(in) :: a,d REAL(8), INTENT(out) :: l,b REAL(8) :: a1,sa,ca,sd,cd REAL(8), PARAMETER :: Leo = 4.936829261 ! 282.85948083 REAL(8), PARAMETER :: L0 = 0.57477039907 ! 32.931918056 REAL(8), PARAMETER :: si = 0.88998807641 ! sin 62.871748611 REAL(8), PARAMETER :: ci = 0.45598379779 ! cos 62.871748611 a1=a-Leo sa=sin(a1); ca=cos(a1) sd=sin(d); cd=cos(d) b=asin(sd*ci-cd*si*sa) l=atan2(sd*si+cd*ci*sa,cd*ca)+L0 END SUBROUTINE Galaxy END MODULE
Листинг E 4.4. Построение распределения звезд по небесной сфере
Program Main USE DFLIB USE HipMain USE Projection IMPLICIT NONE ! Физическое разрешение окна INTEGER, PARAMETER :: MaxX = 1000, MaxY = 500 REAL(8), PARAMETER :: LowX = -2.1 , LowY = -1.05 ! Логические REAL(8), PARAMETER :: HighX= +2.1 , HighY= +1.05 ! координаты REAL(8) :: l, b ! Галактические координаты REAL(8) :: x, y ! Декартовы координаты LOGICAL :: status1 ! Вспомогательные величины INTEGER(2) :: status2 INTEGER(4) :: status4 INTEGER(4) :: color ! Цвет звезды TYPE(THipparcos) :: s TYPE (windowconfig) :: wc ! Свойства графического окна TYPE (wxycoord) :: wxy ! Вспомогательная величина wc.numxpixels = MaxX ! Заполнение структуры свойств окна wc.numypixels = MaxY wc.numtextcols = -1 wc.numtextrows = -1 wc.numcolors = -1 wc.title = "Aitoff"C wc.fontsize = \#0008000C ! 10 X 12 status1 = SETWINDOWCONFIG(wc) ! Инициализация графики IF (.NOT. status1) status1 = SETWINDOWCONFIG(wc) status2=SetWindow(.TRUE.,LowX, LowY, HighX, HighY) status2=INITIALIZEFONTS( ) CALL AitoffGrid(30,.TRUE.) CALL OpenHipparcosMain() ! Открытие каталога и status4=InitCriteria('I-II.txt') ! инициализация критерия do while (ReadHipparcosMain(s)) ! Цикл чтения звезд IF (inCelestia(s.HIP)) THEN ! Проверка критерия SELECT CASE(S.SP(1:1)) ! Определение цвета звезды CASE ('O') color = RGBTOINTEGER(90,64,255) CASE ('B') color = RGBTOINTEGER(128,128,255) CASE ('A') color = RGBTOINTEGER(255,255,255) CASE ('F') color = RGBTOINTEGER(255,255,128) CASE ('G') color = RGBTOINTEGER(255,230,40) CASE ('K') color = RGBTOINTEGER(255,160,0) CASE ('M') color = RGBTOINTEGER(255,0,0) CASE default color=0 END SELECT ! Перевод экваториальных координат в радианы, ! а затем в галактические координаты CALL Galaxy(rad(s.RADeg),rad(s.DEDeg),l,b) ! Вычисление декартовых координат проекции Айтофа CALL Aitoff(l,b,x,y) ! Поставить точку (можно заменить на круг) status4=SetPixelRGB\_w(x,y,color) END IF END DO CALL ClearCriteria() CALL CloseHipparcosMain() ! Сохранить изображение в файле status4=SaveImage("Aitoff.bmp"C,0,0,MaxX-1,Maxy-1) END PROGRAM
Листинг E 4.5. Формирование прямоугольных координат звезд
PROGRAM Main USE HipMain IMPLICIT NONE CHARACTER(*),PARAMETER :: Criteria='O-B5' ! Имя файла критерия INTEGER(4) :: n = 0 ! Счетчик Type(THipparcos) :: s REAL(8) :: r ! Расстояние REAL(8) :: l, b ! Галактические координаты REAL(8) :: x,y,z ! Декартовы галактические координаты INTEGER, PARAMETER :: f = 14 ! Файл вывода результатов OPEN(f, file=Criteria // '.DAT') CALL OpenHipparcosMain() ! Файл списка звезд имеет расширение .TXT PRINT "(I6,' звезд в критерии')",InitCriteria(Criteria//'.txt') do while (ReadHipparcosMain(s)) IF (s.NoPlx) CYCLE ! нет данных о паралл. IF (s.plx<=0.0) CYCLE ! "плохое" значение параллакса IF (s.sigma\_plx/s.plx>0.5) CYCLE ! низкая точность пар. IF (inCelestia(s.HIP)) THEN r=1000.0/s.plx ! Вычисление расстояния в пк IF (r>500.0) CYCLE ! Отброс далеких звезд ! Перевод в галактические координаты CALL Galaxy(rad(s.RADeg),rad(s.DEDeg),l,b) x=r*cos(b)*cos(l) ! Вычисление прямоугольных y=r*cos(b)*sin(l) ! галактических координат z=r*sin(b) write(f,'(3F10.2)'), x,y,z ! Вывод в файл n=n+1 ! Увеличение счетчика на единицу END IF END do CALL ClearCriteria() CALL CloseHipparcosMain() Close(f) print '(I6,"звезд обработано.")',n END PROGRAM
<< 10. Приложение D | Оглавление | Литература >>
Публикации с ключевыми словами:
астрометрия - каталоги - Hipparcos
Публикации со словами: астрометрия - каталоги - Hipparcos | |
См. также:
Все публикации на ту же тему >> |