Algorithms
Alternative Stellar Correlation Theories
A video describing this algorithm is coming soon.
Introduction
The following are two algorithms that, when used together, can generate alternative stellar correlation theories.
Stellarium Script: starwindow.ssc
//======================================================================
// This code is released into the public domain by the author.
//
// You are free to use, modify, distribute,
// and sell this code for any purpose
// without any restrictions or conditions.
// No warranty of any kind is provided.
//
// Author: Andrew Spencer
// Date: June 8, 2023
// File: starwindow.ssc
//======================================================================
//
// Description:
// Script that finds all visible stars in a window surrounding Orion.
// The window of the sky is defined by
// a lower bound on the RA is set by Aldebaran (HIP 21421),
// an upper bound on the RA is set by Sirius (HIP 32349),
// a lower bound on the declenation is -25 degrees, and
// an upper bound on the declenation is +5 degrees.
// The output is saved to file $USERDIR\starwindow.txt
// In Windows $USERDIR is something like
// C:\Users\USERNAME\AppData\Roaming\Stellarium
//
// An example workflow, using Windows with MinGW and Bash's sort, is:
// 1) Run Stellarium script, starwindow.ssc, to generate starwindow.txt.
// (Note: this takes ~10 seconds to run.)
// 2) Move starwindow.txt to GetCandidates.f90's folder.
// 3) gfortran -o GetCandidates.exe -O3 GetCandidates.f90
// 4) GetCandidates.exe
// (Note: this takes ~15 seconds to run.)
// 5) sort -n rawlist.txt > sortedlist.txt
// (Note: this takes less than a second to run.)
// NOTE: Windows sort program doesn't do numerical sorting, which is why
// I use the bash sort.
//
//======================================================================
var maglim=6.5; // upper bound on apparent magnitude of individual stars
// Do not change maglim.
// A more restrictive maxmag can be used in GetCandidates.f90.
core.resetOutput();
core.setObserverLocation(31.13417, 29.97917, 60.0, 1.,
"The Great Pyrmaid", "Earth");
core.setTimeRate(0);
core.setDate("-2574-01-01T00:00:00","local");
core.wait(0.01);
for (var i = 21421; i <= 32349; i++) {
var starMap = core.getObjectInfo("HIP"+i.toString());
var vmag=starMap["vmag"];
if(vmag<=maglim*1.0001) {
var dec=starMap["dec"];
if(dec<=5.0 && dec>=-25.0) {
var ra=starMap["ra"];
core.output(i.toString()+" "+ra.toString()+" "
+dec.toString()+" "+vmag.toString());
}
}
}
core.saveOutputAs("starwindow.txt");
Fortran code: GetCandidates.f90
!=======================================================================
! This code is released into the public domain by the author.
!
! You are free to use, modify, distribute,
! and sell this code for any purpose
! without any restrictions or conditions.
! No warranty of any kind is provided.
!
! Author: Andrew Spencer
! Date: June 8, 2023
! File: GetCandidates.f90
!=======================================================================
!
! This program reads a file generated with Stellarium, starwindow.txt,
! and constructs every possible star triplet and puts them into a
! normalized coordinate system to compare with the Giza Pyramid Complex
! (GPC) layout. For each triplet, error in the correspondance and
! combined magnitude are computed. Finally, the topmost triplets by both
! standards are found and stored in rawlist.txt.
!
! NOTE: This does not sort the candidate star triplets.
! An example workflow, using Windows with MinGW and Bash's sort, is:
! 1) Run Stellarium script, starwindow.ssc, to generate starwindow.txt.
! (Note: this takes ~10 seconds to run.)
! 2) Move starwindow.txt to GetCandidates.f90's folder.
! 3) gfortran -o GetCandidates.exe -O3 GetCandidates.f90
! 4) GetCandidates.exe
! (Note: this takes ~15 seconds to run.)
! 5) sort -n rawlist.txt > sortedlist.txt
! (Note: this takes less than a second to run.)
! NOTE: Windows sort program doesn't do numerical sorting, which is why
! I use the bash sort.
!
! There are three modifications to this source code to consider:
! 1) change the maximum magnitude, maxmag, of every star in triplet,
! 6.5 is "bairly visible" under clearest and darkest skies,
! 3.0 is what I call "bright",
! 2) changing the order in which the errors and magnitudes get printed
! for later sorting (see line 180 where rawlist.txt gets written),
! 3) change the error criteria, r, to be more or less exclusive.
! r of 1/6 is approximately Menkaure's diagonal
! in normalized coordinates
!
! After sorting rawlist.txt, to see the ranking of OCT vs MAAT stars,
! search sortedlist.txt for
! "26727 26311 25930" (the OCT triplet), or
! "27989 25930 24436" (the MAAT triplet)
!=======================================================================
program GetCandidates
implicit none
double precision, parameter :: maxmag=6.500001d0 !3.000001d0
double precision, parameter :: r=1.d0/6.d0
integer :: nstars
integer, dimension(:), allocatable :: hip_index
double precision, dimension(:), allocatable::ra_arr,dec_arr,vmag,pmag
integer :: triplets
!hips store the HIP indices of candidate star triplets
integer, dimension(:,:), allocatable :: hips
double precision, dimension(:), allocatable :: errs,cmags,mmags
double precision :: ra,dec,x1,y1,x2,y2,x3,y3,norm,mag1,mag2,mag3
double precision :: coords(3,3),com(3),sc1(3),sc2(3)
double precision :: ang,cang,sang,Menx,Meny,dx,dx2,dy,dy2,r2
double precision :: dyt,dyt2,dist2,x3t,y3t
double precision, parameter :: degtorad = acos(-1.d0)/180.d0
integer :: i, j, k, n, status
call get_line_count("starwindow.txt", nstars)
triplets=nstars*(nstars-1)*(nstars-2)
allocate(hip_index(nstars))
allocate(ra_arr(nstars),dec_arr(nstars),vmag(nstars),pmag(nstars), &
hips(3,triplets),errs(triplets),cmags(triplets),mmags(triplets))
hips=0;errs=0.d0;cmags=0.d0;mmags=0.d0
open(unit=10,file="starwindow.txt",status='old',action='read')
do i = 1, nstars
read(10, *) hip_index(i), ra_arr(i), dec_arr(i), vmag(i)
pmag(i)=10.d0**(-0.4d0*vmag(i))
ra_arr(i)=ra_arr(i)*degtorad
dec_arr(i)=dec_arr(i)*degtorad
enddo
close(10)
!Get norm and normalized coordinates of Menkaure in GPC, (Menx,Meny)
x1=31.13417; y1=29.97917 ! Khufu's pyramid
x2=31.13083; y2=29.97611 ! Khafre's pyramid
x3=31.12833; y3=29.9725 ! Menkaure's pyramid
x2=x2-x1; y2=y2-y1 !Put Khufu at (0,0)
norm = dsqrt(x2*x2+y2*y2)!Distance between Khufu and Khafre is norm
x2=x2/norm; y2=y2/norm !Change Khafre to a unit vector
ang=datan2(y2,x2) !Angle from x-axis to Khafre's unit vector
cang=dcos(ang); sang=dsin(ang)
x3=(x3-x1)/norm; y3=(y3-y1)/norm !Normalize Menkaure
Menx= x3*cang+y3*sang !Rotate Menkaure so Khafre is at (0,1)
Meny=-x3*sang+y3*cang !Rotate Menkaure so Khafre is at (0,1)
!Set the error criterion
r2=r*r
n=0 !n is the count of stars meeting the criteria
do i = 1, nstars
mag1 = vmag(i)
if(mag1>maxmag) cycle
ra = ra_arr(i); dec = dec_arr(i)
coords(1,1) = dcos(dec)*dcos(ra)
coords(2,1) = dcos(dec)*dsin(ra)
coords(3,1) = dsin(dec)
do j = 1, nstars
mag2 = vmag(j)
if(j==i.or.mag2>maxmag) cycle
ra = ra_arr(j); dec = dec_arr(j)
coords(1,2) = dcos(dec)*dcos(ra)
coords(2,2) = dcos(dec)*dsin(ra)
coords(3,2) = dsin(dec)
do k = 1, nstars
mag3 = vmag(k)
if(k==i.or.k==j.or.mag3>maxmag) cycle
ra = ra_arr(k); dec = dec_arr(k)
coords(1,3) = dcos(dec)*dcos(ra)
coords(2,3) = dcos(dec)*dsin(ra)
coords(3,3) = dsin(dec)
!Rotate the triplet so the center of mass points to (0,0,-1),
!then project to 2D using the stereographic projection
com = sum(coords,2)
com = com/dsqrt(sum(com*com))
sc1(1)= com(2)*com(2)/(1.d0-com(3)) - com(3)
sc1(2)=-com(1)*com(2)/(1.d0-com(3))
sc1(3)=com(1)
sc2(1)=sc1(2)
sc2(2)= com(1)*com(1)/(1.d0-com(3)) - com(3)
sc2(3)=com(2)
norm=1.d0+sum(com*coords(:,1))
x1 = sum(sc1*coords(:,1))/norm
y1 = sum(sc2*coords(:,1))/norm
norm=1.d0+sum(com*coords(:,2))
x2 = sum(sc1*coords(:,2))/norm
y2 = sum(sc2*coords(:,2))/norm
norm=1.d0+sum(com*coords(:,3))
x3 = sum(sc1*coords(:,3))/norm
y3 = sum(sc2*coords(:,3))/norm
!After projection is finished,
!put triplet into normlized coordinates.
!Repeat coordinate transform for the GPC,
!but now for star triplets,
x2=x2-x1; y2=y2-y1
norm = dsqrt(x2*x2+y2*y2)
x2=x2/norm; y2=y2/norm
ang=datan2(y2,x2)
cang=dcos(ang); sang=dsin(ang)
x3t=(x3-x1)/norm; y3t=(y3-y1)/norm
x3= x3t*cang+y3t*sang
y3=-x3t*sang+y3t*cang
!Find candidates for Giza triplet
dx=x3-Menx; dy=y3-Meny; dyt=-y3-Meny
dx2=dx*dx; dy2=dy*dy; dyt2=dyt*dyt
dist2=min(dx2+dy2,dx2+dyt2)
if(dist2<=r2) then !Check if triplet meets criteria
!If found, store info in arrays
n=n+1 !Keep count of triplets meeting the criteria
errs(n)=dsqrt(dist2)
cmags(n)=-2.5d0*dlog10(pmag(i)+pmag(j)+pmag(k))
mmags(n)=max(mag1,mag2,mag3)
hips(1,n)=hip_index(i)
hips(2,n)=hip_index(j)
hips(3,n)=hip_index(k)
endif
enddo
enddo
enddo
deallocate(hip_index,ra_arr,dec_arr,vmag,pmag)
open(unit=10,file="rawlist.txt",status='unknown', &
action='write',iostat=status)
do i=1,n !put errors or combined magnitudes first for later sorting
write(10, '(3(F15.10,1X), 3(I5,1X))') &
errs(i),cmags(i),mmags(i),hips(:,i)
enddo
close(10)
deallocate(hips,errs,cmags,mmags)
contains
!=======================================================================
! get_line_count serves two functions:
! 1) counts the number of lines in a text file, and
! 2) stops execution at beginning of run if problem reading file.
!=======================================================================
subroutine get_line_count(filename, line_count)
implicit none
character(len=*), intent(in) :: filename
integer, intent(out) :: line_count
integer :: iostat
character(len=256) :: line
! Initialize variables
line_count = 0
! Open the file
open(unit=11,file=filename,status='old',action='read',iostat=iostat)
if (iostat /= 0) then
write(*,*) "Error opening the file."; stop
endif
! Read the file line by line and count the lines
do
read(11, '(A)', iostat=iostat) line
if (iostat /= 0) return
line_count = line_count + 1
end do
! Close the file
close(11)
end subroutine get_line_count
end program GetCandidates
Comments
Post a Comment