APSF
//- Variable Star Ephemeris. This script works out the dates/times of the minima of variable star objects.
//- It produces a tabular report for a given period.
//-
//- Paul Rodman, June 2008 (with input from David O'Driscoll)
//-
//- NOTES:
//- a. Only objects in the plan that originated from the GCVS catalogue will be considered, 
//-    since their Notes field contains required information to compute the times of minima.
//- b. The script assumes you have set up two Site resources to reflect the correct location
//-    of your Home and your Observatory (these can be the same site if you observe from home,
//-    or different if you observe from a remote Observatory site.
//- c. The script assumes that you have set your plan document to the Observatory site before 
//-    starting (since there is currently no method of changing that from within a script).
//- 
//- Version 1.0  26 June 2008
//-   Initial release (Note: does not take into account visibility at Observatory site)
//- Version 1.1  26 June 2008
//-   Fixed bug with "Month from specified date" option
//-   Fixed bug with malformed Notes fields
//-   Added visibility option
//-   Added ability to show object altitude
//- Version 1.2  27 June 2008
//-   Fixed bug with site "sanity check"
//- Version 1.3  4 July 2008
//-   Added table orientation feature to flip rows and columns.

//------------------------------------------------------------------------------------------
function Normalise(t as double, mint as double, maxt as double) as double
  dim r as double
  r=t
  while r<mint
    r=r+(maxt-mint)
  wend
  while r>=maxt
    r=r-(maxt-mint)
  wend
  return r
end function

//------------------------------------------------------------------------------------------
function Julian(year as integer, month as integer, day as double) as double
  // Meeus Pg 60-61
  
  dim y,m,a,b as integer
  dim d as double
  dim JulianCalendar as boolean
  
  y=year
  m=month
  d=day
  if month<=2 then
    y=y-1
    m=m+12
  end if
  JulianCalendar = year<1582 or (year=1582 and ((month<10) or month=10 and day<5))
  if JulianCalendar then
    b=0
  else
    a=floor(y/100)
    b=2-a+floor(a/4)
  end if
  return floor(365.25*(y+4716))+floor(30.6001*(m+1))+d+b-1524.5
end function

//------------------------------------------------------------------------------------------
function Julian(dt as double) as double
  return Julian(YearOfDate(dt),MonthOfDate(dt),DayOfDate(dt)+HourOfDate(dt)/24.0+MinuteOfDate(dt)/1440.0+SecondOfDate(dt)/86400.0)
end function

//------------------------------------------------------------------------------------------
function GetMinima(ut as double, epoch as double, period as double) as double()
  dim n as integer, jd,minima(-1),t as double
  jd=Julian(ut)
  n=ceil((jd-epoch)/period)
  t=epoch+n*period
  for i=0 to 10000
    if t>=jd+1.0 then exit
    minima.Append ut+(t-jd)*86400.0
    t=t+period
  next
  return minima
end function

//------------------------------------------------------------------------------------------
function DoubleToDateTime(d as double) as string
  return DoubleToDate(d)+" "+DoubleToTime(d,false)
end function

//------------------------------------------------------------------------------------------
function TimeInRange(dt as double, tfrom as double, tto as double) as boolean
  dim t as double
  t=HourOfDate(dt)+MinuteOfDate(dt)/60.0+SecondOfDate(dt)/3600.0
  if tfrom<tto then
    return t>=tfrom and t<=tto
  else
    return t>=tfrom or t<=tto
  end if
end function

//------------------------------------------------------------------------------------------
function IsVisible(ob as APPlanObject, t as double, minalt as double) as boolean
  PlanLocalDateTime=t
  return ob.Altitude>=minalt
end function

//------------------------------------------------------------------------------------------
function ExtractData(ob as APPlanObject, byref epoch as double, byref period as double) as boolean
   
  dim j,k,l as integer, s,c as string
   
  s=ob.Notes
  j=Instr(s,"Epoch:")
  if j<=0 then return false
  j=j+6
  l=len(s)
  for k=j to len(s)
    c=mid(s,k,1)
    if not IsNumeric(c) and c<>"." then
      l=k-1
      exit
    end if
  next
  epoch=val(mid(s,j,l-j+1))
  if epoch<2.4E6 then return false
  
  j=Instr(s,"Period:")
  if j<=0 then return false
  j=j+7
  l=len(s)
  for k=j to len(s)
    c=mid(s,k,1)
    if not IsNumeric(c) and c<>"." then
      l=k-1
      exit
    end if
  next
  period=val(mid(s,j,l-j+1))
  if period<=0.0 then return false
  return true  
end function

//------------------------------------------------------------------------------------------

dim s,ss(-1),srow(-1),fontlist(-1),stitle as string
dim timeperiod,which,homesite,obssite,showtime(2),twrt,i,j,k,d,m,y,ndays,iob,row,col as integer
dim startdate,tfrom,tto,minalt,epoch(-1),period(-1),diff,dt,diffObstoUT,diffHometoUT,ep,pe,minima(-1),t,rowhgt,rhgt(-1),savet as double
dim ob,oblist(-1) as APPlanObject
dim ok,printreport,doVisibility,showAltitude as boolean
dim c as Canvas, tbl,tblhdr as Table
dim pfont,dfont,fontlist(-1) as string, psize,dsize,orientation,startday,endday as integer
dim tstart as double

const maxcols = 5
const maxdaycols = 11

fontlist=Fonts
if ubound(fontlist)<0 then
  print "No fonts defined in your OS!"
  return
end if

// Load persistent values
SaveRestoreGlobal(true)

timeperiod=RestoreIntegerValue("timeperiod",0)
which=RestoreIntegerValue("which",0)
homesite=RestoreIntegerValue("homesite",0)
obssite=RestoreIntegerValue("obssite",0)
showtime(0)=1
showtime=RestoreIntegerValue("showtime",showtime)
tfrom=RestoreDoubleValue("tfrom",20.0)
tto=RestoreDoubleValue("tto",5.0)
twrt=RestoreIntegerValue("twrt",1)
minalt=RestoreDoubleValue("minalt",5.0)
pfont="Times"
if fontlist.IndexOf(pfont)<0 then 
  pfont="Times New Roman"
end if
if fontlist.IndexOf(pfont)<0 then 
  pfont="System"
end if
if fontlist.IndexOf(pfont)<0 then 
  pfont=fontlist(0)
end if
dfont="System"
if fontlist.IndexOf(dfont)<0 then 
  dfont=fontlist(0)
end if
pfont=RestoreStringValue("pfont",pfont)
dfont=RestoreStringValue("dfont",dfont)
psize=RestoreIntegerValue("psize",10)
dsize=RestoreIntegerValue("dsize",11)
doVisibility=RestoreBooleanValue("doVisibility",true)
showAltitude=RestoreBooleanValue("showAltitude",true)
orientation=RestoreIntegerValue("orientation",0)

startdate=CurrentDate

// Set up dialog

redim ss(-1)
ss.Append "Day"
ss.Append "Week"
ss.Append "Month (from 1st day)"
ss.Append "Month (from specified day)"
SetPopupParameter("Time Period",timeperiod,ss)
SetDateParameter(true,"Start Date (at Home site)",startdate)
SetChoiceParameter("Objects to consider",which,"All objects","Selected object","Highlighted objects")
SetChoiceParameter(true,"Table orientation",orientation,"Rows: Date, Columns: Objects", "Rows: Objects, Columns: Date")
SetCheckListParameter("Show time of minima for",showtime,"Home site","Observatory site","UT")
SetBooleanParameter(true,"Show object altitude",showAltitude)
SetPopupParameter("Home Site",homesite,Sites)
SetPopupParameter(true,"Observatory Site",obssite,Sites)
SetTimeParameter("Observing session time From",tfrom)
SetTimeParameter(true,"Observing session time To",tto)
SetChoiceParameter("Observing session time with respect to",twrt,"Home site","Observatory site","UT")
SetBooleanParameter("Consider object visibility",doVisibility)
SetDoubleParameter(true,"Minimum altitude (deg)",minalt,0.0,90.0)
SetPopupParameter("Display Font",fontlist.IndexOf(dfont),fontlist)
SetIntegerParameter(true,"Size.1",dsize,4,128)
SetPopupParameter("Print Font",fontlist.IndexOf(pfont),fontlist)
SetIntegerParameter(true,"Size.2",psize,4,128)

ParameterWindowOKCaption("Display")
AddPushbutton("Print...")

if not EditParameters("Variable Star Ephemeris") then return

// Retrieve dialog values
timeperiod=GetPopupParameter("Time Period")
startdate=GetDateParameter("Start Date (at Home site)")
which=GetChoiceParameter("Objects to consider")
orientation=GetChoiceParameter("Table orientation")
homesite=GetPopupParameter("Home Site")
obssite=GetPopupParameter("Observatory Site")
showtime=GetCheckListParameter("Show time of minima for")
showAltitude=GetBooleanParameter("Show object altitude")
tfrom=GetTimeParameter("Observing session time From")
tto=GetTimeParameter("Observing session time To")
twrt=GetChoiceParameter("Observing session time with respect to")
doVisibility=GetBooleanParameter("Consider object visibility")
minalt=GetDoubleParameter("Minimum altitude (deg)")
dfont=fontlist(GetPopupParameter("Display Font"))
dsize=GetIntegerParameter("Size.1")
pfont=fontlist(GetPopupParameter("Print Font"))
psize=GetIntegerParameter("Size.2")

printreport=GetPushButton("Print...")

// Save persistent values

SaveIntegerValue("timeperiod",timeperiod)
SaveIntegerValue("which",which)
SaveIntegerValue("homesite",homesite)
SaveIntegerValue("obssite",obssite)
SaveIntegerValue("showtime",showtime)
SaveDoubleValue("tfrom",tfrom)
SaveDoubleValue("tto",tto)
SaveIntegerValue("twrt",twrt)
SaveDoubleValue("minalt",minalt)
SaveStringValue("pfont",pfont)
SaveStringValue("dfont",dfont)
SaveIntegerValue("psize",psize)
SaveIntegerValue("dsize",dsize)
SaveBooleanValue("doVisibility",doVisibility)
SaveBooleanValue("showAltitude",showAltitude)
SaveIntegerValue("orientation",orientation)

if FixedDate then 
  savet = PlanLocalDateTime
  PlanLocalDateTime = 0
else 
  savet = 0
end if

// Do site setting sanity check
dt = (PlanLocalDateTime-Site(obssite+1).GMTOffset(PlanLocalDateTime)*3600.0)
diff=abs(HourOfDate(dt)+MinuteOfDate(dt)/60.0+SecondOfDate(dt)/3600.0 - PlanGMT)*60.0
if diff>=1.0 then
  PlanLocalDateTime=savet
  print "Current site does not appear to be set to the observatory site: "+Site(obssite+1).Name
  return
end if

// Find usable objects
for i=1 to nObj
  ob=Obj(i)
  select case which
  case 0 // All objects
    ok=true
  case 1 // Selected object
    ok=ob.Selected
  case 2 // Highlighted
    ok=ob.IsHighlighted
  end select
  if ok then
    // Check and see if Notes contains the requisite information
    if ExtractData(ob,ep,pe) then
      oblist.Append ob
      epoch.Append ep
      period.Append pe
    end if
  end if
next

if ubound(oblist)<0 then
  PlanLocalDateTime=savet
  print "No applicable variable star objects found!"
  return
end if

// Figure out date range
select case timeperiod
case 0 // Day
  ndays=1
case 1 // Week
  ndays=7
case 2 // Month, from 1st
  d=1
  m=MonthOfDate(startdate)
  y=YearOfDate(startdate)
  ndays=DaysInMonth(m,y)
  startdate=MakeDate(y,m,d)
case 3 // Month, from specified
  d=DayOfDate(startdate)
  m=MonthOfDate(startdate)
  y=YearOfDate(startdate)
  ndays=DaysInMonth(m,y)-d+1
  if d>1 then
    m=m+1
    if m>12 then
      m=1
      y=y+1
    end if
    ndays=ndays+min((d-1),DaysInMonth(m,y))
  end if
end select

// Figure out session times at observatory
// Assume DST same as that for first day

diffObstoUT=-Site(obssite+1).GMTOffset(startdate)
diffHometoUT=-Site(homesite+1).GMTOffset(startdate)

select case twrt
case 0 // Home site
  tfrom=Normalise(tfrom+diffHomeToUT-diffObsToUT,0.0,24.0)
  tto=Normalise(tto+diffHomeToUT-diffObsToUT,0.0,24.0)
case 1 // Observatory site
  // No change required
case 2 // UT
  tfrom=Normalise(tfrom-diffObsToUT,0.0,24.0)
  tto=Normalise(tto-diffObsToUT,0.0,24.0) 
end select

// Set up canvas
stitle="Variable Star Planner for "+DoubleToDate(startdate)+" to "+DoubleToDate(startdate+(ndays-1)*86400.0)

if printreport then
  c=new Canvas(false)
  c.TextFont(pfont,psize)
else
  c=new Canvas(1100,800,stitle)
  c.TextFont(dfont,dsize)
end if

if c.Cancelled then
  c.Close
  return
end if

if orientation=0 then
  // Loop for each object
  col=1
  for iob=0 to ubound(oblist) 
    ob=oblist(iob)
  
    if tbl=nil or col>maxcols then
      if tbl<>nil then 
        y=c.DrawTable(tbl,0,0,c.Width,c.Height,grid_Thin)
        c.NewPage
      end if
      tbl=new Table(ndays,maxcols+1,stitle)
      tbl.ColumnWidth(1)=7
      tbl.ColumnTitle(1)="Date"+EndOfLine+"(Home)"
      tbl.RowHeight(0)=200.0
  	   tbl.RowStyle(0)=style_Bold+style_Inverted+Style_Gray
      t=startdate
      redim rhgt(ndays)
      for d=1 to ndays
        rhgt(d)=100.0
        tbl.Cell(d,1)=DoubleToDate(t)
        t=t+86400.0
      next
      col=1
    end if
  
    s=ob.ID
    if ob.Name<>"" and ob.Name<>ob.ID then s=s+EndofLine+ob.Name
    tbl.ColumnTitle(col+1)=s
    
    // Loop through the days
    t=startdate
    for d=1 to ndays
      minima=GetMinima(t+diffHometoUT*3600.0,epoch(iob),period(iob))
      redim srow(-1)
      for i=0 to ubound(minima)
        if TimeInRange(minima(i)-diffObstoUT*3600.0,tfrom,tto) then
          if not doVisibility or IsVisible(ob,minima(i)-diffObstoUT*3600.0,minalt) then
            redim ss(-1)
            if showtime(0)<>0 then ss.Append "Home: "+DoubleToTime(minima(i)-diffHometoUT*3600.0,false)
            if showtime(1)<>0 then ss.Append "Obs: "+DoubleToTime(minima(i)-diffObstoUT*3600.0,false)
            if showtime(2)<>0 then ss.Append "UT: "+DoubleToTime(minima(i),false)
            if showAltitude then ss.Append "Alt: "+str(round(ob.Altitude))+DegreeSymbol
            s=Join(ss,", ")
            if ubound(ss)=0 then s=NthField(s,": ",2)
            srow.Append s
          end if
        end if
      next
      rhgt(d)=max(rhgt(d),(ubound(srow)+1)*100.0)
      tbl.RowHeight(d)=rhgt(d)
      tbl.Cell(d,col+1)=Join(srow,EndOfLine)
      t=t+86400.0
    next
    col=col+1
  next

  y=c.DrawTable(tbl,0,0,c.Width,c.Height,grid_Thin)
else
  startday=1
  endday=ndays
  tstart=startdate
  while startday<=endday
  if tstart>startdate then c.NewPage
  ndays=min(maxdaycols,endday-startday+1)
  // Create header and data tables
  tblhdr=new Table(1,ndays+1,stitle)
  tbl=new Table(1,ndays+1)
  tblhdr.ColumnWidth(1)=10
  tbl.ColumnWidth(1)=10
  tblhdr.ColumnTitle(1)="Star"
  tbl.RowHeight(1)=200.0
  t=tstart
  for i=1 to ndays
    tblhdr.ColumnTitle(i+1)=DoubleToDate(t)
    t=t+86400.0
  next
  tblhdr.RowStyle(0)=style_Bold+style_Inverted+Style_Gray
  tblhdr.DeleteRow(1)
  y=0
  // Loop for each object
  for iob=0 to ubound(oblist) 
    ob=oblist(iob)
    s=ob.ID
    if ob.Name<>"" and ob.Name<>ob.ID then s=s+EndofLine+ob.Name
    tbl.Cell(1,1)=s
    
    // Loop through the days
    rowhgt=200.0
    t=tstart
    for d=1 to ndays
      minima=GetMinima(t+diffHometoUT*3600.0,epoch(iob),period(iob))
      redim srow(-1)
      for i=0 to ubound(minima)
        if TimeInRange(minima(i)-diffObstoUT*3600.0,tfrom,tto) then
          if not doVisibility or IsVisible(ob,minima(i)-diffObstoUT*3600.0,minalt) then
            redim ss(-1)
            if showtime(0)<>0 then ss.Append "Home: "+DoubleToTime(minima(i)-diffHometoUT*3600.0,false)
            if showtime(1)<>0 then ss.Append "Obs: "+DoubleToTime(minima(i)-diffObstoUT*3600.0,false)
            if showtime(2)<>0 then ss.Append "UT: "+DoubleToTime(minima(i),false)
            if showAltitude then ss.Append "Alt: "+str(round(ob.Altitude))+DegreeSymbol
            s=Join(ss,EndOfLine)
            if ubound(ss)=0 then s=NthField(s,": ",2)
            srow.Append s
          end if
        end if
      next
      s=Join(srow,EndOfLine)
      rowhgt=max(rowhgt,CountFields(s,EndOfLine)*100.0)
      tbl.RowHeight(1)=rowhgt
      tbl.Cell(1,d+1)=s
      t=t+86400.0
    next
    
    if y>c.Height*0.9 then 
      y=0
      c.NewPage
    end if
    if y=0 then y=c.DrawTable(tblhdr,0,0,c.Width,c.Height,grid_Thin)
    y=c.DrawTable(tbl,0,y,c.Width,c.Height,grid_Thin)
  next
  startday=startday+maxdaycols
  tstart=tstart+maxdaycols*86400.0
  wend
end if

c.Close

PlanLocalDateTime=savet