APSF
//- Find guide stars close to the selected object, and add them to the plan.
//- 
//- Paul Rodman, Jan 2008
//- Bruce Pipes (idea and debugging)
//- 
//- V1.0  22 Jan 2008
//-   Initial release
//- V1.1  22 Jan 2008
//-   Ability to select guide stars to add to plan
//- V1.2  23 Jan 2008
//-   Corrected a sorting problem (note: will require V1.6.1 to function correctly)
//- V1.3  24 Jan 2008
//-   Added bearing info for guide stars
//- V1.3.1  24 Jan 2008
//-   Fixed problem with bearing in notes field
//- V1.4  24 Jan 2008
//-   Additional option to choose direction of bearing wrt N
//- V1.5  4 Feb 2008
//-   Make the catalogue selector a list, and add magnitude info
//- V1.6  12 Feb 2008  (Requires application V1.6.1)
//-   Allow for reticle/telescope/aid info to be used to compute search annulus

//# Author: Paul Rodman
//# E-Mail: paul@ilanga.com
//# Name: Guide Stars
//# Category: Imaging|Guide Stars
//# MinVersion: 1.6.1
//# Keywords: Imaging, Guide Stars

//----------------------------------------------------------------------------------------
Function BearingTo(ra1 as double, dec1 as double, ra2 as double, dec2 as double) As double
  dim a,g,b,c,alpha,diff,p as double
  
  const pi=3.141592654
  const DegreesToRadians = 0.017453292
  const HoursToRadians = 0.261799387
  const RadiansToDegrees = 57.29577951
  
  if ra1=ra2 and dec1=dec2 then return 0.0
  
  b=(90.0-dec1)*DegreesToRadians
  c=(90.0-dec2)*DegreesToRadians
  alpha=abs(ra1-ra2)*HoursToRadians
  if alpha>pi then
    alpha=2.0*pi-alpha
  end if
  g=acos(cos(b)*cos(c)+sin(b)*sin(c)*cos(alpha))
  
  p=cos(dec1*DegreesToRadians)*sin((ra1-ra2)*HoursToRadians)/sin(g)
  if p<=-1 then return 90.0
  if p>=1 then return 270.0
  
  a=asin(p)
  
  if dec2<dec1 then
    a=a+pi
  else
    if ra2-ra1<=-12.0 or ra2>ra1 then a=-a else a=2.0*pi-a
  end if
  
  a=a*RadiansToDegrees
  if a<0.0 then a=a+360.0
  if a>=360.0 then a=a-360.0
  return a
  
End Function

//----------------------------------------------------------------------------------------
sub ComputeDistances(r as APReticle, t as TelescopeResource, a as VisualAidResource, byref distFrom as double, byref distTo as double)
  dim flength as double
  
  const RadiansToDegrees = 57.29577951

  if a=nil then
    flength=t.FocalLength
  else
    flength=(t.FocalLength+a.FocalLength)*a.Magnification
  end if

  distFrom=max(r.Offset-0.5*r.Height,0.0)/flength*RadiansToDegrees
  distTo=(r.Offset+0.5*r.Height)/flength*RadiansToDegrees
end sub

//----------------------------------------------------------------------------------------
dim srch as APSearch, target,ob as APPlanObject, i,j,k,btype as integer, gstars(-1) as APCatalogObject, cats(-1) as APCatalog
dim t,t1,t2,mag as double, sz,sel(-1),u4(-1),targetID,catname(-1) as string, selb(-1),where as integer, bearing as boolean
dim searchCat as integer, magFrom, magTo, distFrom, distTo as double, targetRA, targetDec, mfrom, mto as double
dim ret,scope,aid as integer, retname(-1),aids(-1) as string, ar as VisualAidResource, rets(-1) as APReticle

// Check than an object is selected
if SelectedObject<1 then
  print "No object selected in plan!"
  return
end if
target=Obj(SelectedObject)
targetRA=target.RA
targetDec=target.Dec
targetID=target.ID

SaveRestoreGlobal(true)
SaveRestoreTag("GuideStars")

// Restore values saved from last run
searchCat = RestoreIntegerValue("searchCat",0)
ret = RestoreIntegerValue("reticle",0)
scope = RestoreIntegerValue("scope",0)
aid = RestoreIntegerValue("aid",0)
btype = RestoreIntegerValue("btype",0)
magFrom = RestoreDoubleValue("magFrom",2.0)
magTo = RestoreDoubleValue("magTo",10.0)
distFrom = RestoreDoubleValue("distFrom",0.25)
distTo = RestoreDoubleValue("distTo",0.5)

// Use dialog to setup parameter values
cats=CatalogsWithType(type_Star)
for i=0 to ubound(cats)
  cats(i).MagnitudeRangeOfType(type_Star,mfrom,mto)
  catname.Append cats(i).Name+"   ("+format(cats(i).nObjectsOfType(type_Star),",0")+" stars, Mag "+ _
    format(mto,"-0.0")+" to "+format(mfrom,"-0.0")+")"
next
SetListChoiceParameter("Search catalogue",searchCat,catname,false,10)
SetDoubleParameter("Magnitude From",magFrom,0.0,25.0)
SetDoubleParameter("Magnitude To",magTo,0.0,25.0)
if target.Size="" then
  sz="Target ("+target.ID+") size is unknown"
else
  sz=trim(NthField(target.Size,"x",1))
  select case right(sz,1)
  case """"
    t1=val(left(sz,len(sz)-1))/3600.0
  case "'"
    t1=val(left(sz,len(sz)-1))/60.0
  case DegreeSymbol,"d"
    t1=val(left(sz,len(sz)-1))
  else
    t1=val(sz)/60.0
  end select
  sz=trim(NthField(target.Size,"x",2))
  if sz="" then
    t2=0.0
  else
    select case right(sz,1)
    case """"
      t2=val(left(sz,len(sz)-1))/3600.0
    case "'"
      t2=val(left(sz,len(sz)-1))/60.0
    case DegreeSymbol,"d"
      t2=val(left(sz,len(sz)-1))
    else
      t2=val(sz)/60.0
    end select
  end if  
  sz=format(t1,"0.00")+DegreeSymbol
  if t2>0.0 then sz=sz+" x "+format(t2,"0.00")+DegreeSymbol
      
  sz="Target ("+target.ID+") size is "+sz+" (Radius: "+format(max(t1,t2)*0.5,"0.00")+DegreeSymbol+")"
end if
SetCaptionParameter(sz,1,true)
retname.Append "Use distances below"
for i=1 to ReticleCount
  if Reticle(i).Shape=2 then 
    rets.Append Reticle(i)
    retname.Append Reticle(i).Name
  end if
next
SetPopupParameter("Reticle",ret,retname)
SetPopupParameter("Telescope",scope,Telescopes)
aids.Append "None"
for i=1 to nVisualAids
  aids.Append VisualAid(i).Name
next
SetPopupParameter("Visual Aid",aid,aids)

SetDoubleParameter("Minimum distance from target "+DegreeSymbol,distFrom,0.0,5.0)
SetDoubleParameter("Maximum distance from target "+DegreeSymbol,distTo,0.0,15.0)
SetChoiceParameter("Bearings",btype,"Clockwise (From N through W)","Counter-clockwise (From N through E)")

ParameterDependency("Minimum distance from target "+DegreeSymbol,"-Reticle")
ParameterDependency("Maximum distance from target "+DegreeSymbol,"-Reticle")
ParameterDependency("Telescope","Reticle")
ParameterDependency("Visual Aid","Reticle")

if not EditParameters then return

searchCat=GetListChoiceParameter("Search catalogue")
magFrom=GetDoubleParameter("Magnitude From")
magTo=GetDoubleParameter("Magnitude To")
ret=GetPopupParameter("Reticle")
scope=GetPopupParameter("Telescope")
aid=GetPopupParameter("Visual Aid")
distFrom=GetDoubleParameter("Minimum distance from target "+DegreeSymbol)
distTo=GetDoubleParameter("Maximum distance from target "+DegreeSymbol)
btype=GetChoiceParameter("Bearings")

if magFrom>magTo then
  t=magFrom
  magFrom=magTo
  magTo=t
end if

if distFrom>distTo then
  t=distFrom
  distFrom=distTo
  distTo=t
end if

// Save parameter values for next run
SaveIntegerValue("searchCat",searchCat)
SaveIntegerValue("reticle",ret)
SaveIntegerValue("scope",scope)
SaveIntegerValue("aid",aid)
SaveIntegerValue("btype",btype)
SaveDoubleValue("magFrom",magFrom)
SaveDoubleValue("magTo",magTo)
SaveDoubleValue("distFrom",distFrom)
SaveDoubleValue("distTo",distTo)

// Compute distFrom and distTo if reticle used

if ret>0 then 
  if aid>0 then ar=VisualAid(aid) else ar=nil
  ComputeDistances(rets(ret-1),Telescope(scope+1),ar,distFrom,distTo)
end if

// Set up catalogue search
srch = new APSearch

srch.AddBound(item_Magnitude, rel_GE, magFrom)
srch.AddBound(item_Magnitude, rel_LE, magTo)
srch.AddType(type_Star)
srch.SetCatalog(cats(searchCat).Name)
srch.WithinDistance(distTo,target.RA,target.Dec)

// Perform search
if not srch.Search then
  print "Error: "+srch.Error
  return
end if

// Extract useful stars from the results
for i=1 to srch.nObjects
  if AngleBetween(target, srch.GetObject(i))>=distFrom then gstars.Append srch.GetObject(i)
next

if ubound(gstars)<0 then
  print "No suitable stars found!"
  return
end if

mag=999.0
j=-1
for i=0 to ubound(gstars)
  t=BearingTo(target.RA,target.Dec,gstars(i).RA,gstars(i).Dec)
  if btype=0 then t=360.0-t
  sel.Append gstars(i).ID+"  Mag: "+format(gstars(i).Magnitude,"-0.0")+" Bearing: "+ _
  format(t,"-0.0")+DegreeSymbol
  if gstars(i).Magnitude<mag then
    j=i
    mag=gstars(i).Magnitude
  end if
next
redim selb(ubound(sel))
selb(j)=1

where = RestoreIntegerValue("where",0)
bearing = RestoreBooleanValue("bearing",false)

sz="Search annulus: "+format(distFrom,"-0.00")+DegreeSymbol+" to "+format(distTo,"-0.00")+DegreeSymbol
SetCaptionParameter(sz,1,true)
SetCheckListParameter("Select stars to add to plan",selb,sel,false,15)
SetChoiceParameter("Where to add the stars",where,"Before the selected object","After the selected object","At the end of the plan")
SetBooleanParameter("Add bearing info to Notes",bearing)

if not EditParameters("Selected Stars ("+str(ubound(sel)+1)+")") then return

selb=GetCheckListParameter("Select stars to add to plan")
where=GetChoiceParameter("Where to add the stars")
bearing=GetBooleanParameter("Add bearing info to Notes")

SaveIntegerValue("where",where)
SaveBooleanValue("bearing",bearing)

// Remove unwanted stars
for i=ubound(gstars) downto 0
  if selb(i)=0 then gstars.Remove i
next

if ubound(gstars)<0 then
  print "No stars selected!"
  return
end if

// record current value of User(4)
redim sel(nObj-1)
redim u4(nObj-1)
for i=1 to nObj
  sel(i-1)=Obj(i).ID+"|"+Obj(i).Name
  u4(i-1)=Obj(i).User(4)
next

// Number the objects
j=1
for i=1 to nObj
  if i=SelectedObject then
    select case where
    case 0  // Before
      k=j
      j=j+ubound(gstars)+1
      Obj(i).User(4)=format(j,"00000")
    case 1 // After
      Obj(i).User(4)=format(j,"00000")
      k=j+1
      j=j+ubound(gstars)+1
    else  // End of plan
      Obj(i).User(4)=format(j,"00000")
      k=nObj+1
    end select
  else
    Obj(i).User(4)=format(j,"00000")
  end if
  j=j+1
next
  
// Add the stars to the plan
for i=0 to ubound(gstars)
  ob=gstars(i).AddToPlan
  ob.User(4)=format(k,"00000")
  if bearing then 
    t=BearingTo(targetRA,targetDec,gstars(i).RA,gstars(i).Dec)
    if btype=0 then t=360.0-t
    ob.Notes=ob.Notes+"[Bearing w.r.t. "+targetID+": "+ _
      format(t,"-0.0")+DegreeSymbol+"]"
  end if
  k=k+1
next

SortObservingList(list_User4)
SortObservingList()

// Recover old User(4) values
for i=1 to nObj
  j=sel.IndexOf(Obj(i).ID+"|"+Obj(i).Name)
  if j>=0 then Obj(i).User(4)=u4(j) else Obj(i).User(4)=""
next