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