APSF
//- Find the opposition point (the point in the sky opposite the sun along the ecliptic).
//- Add it as an object to the plan
//-
//- Paul Rodman, Jan 2010
//-
//- V1.0   27 Jan 2010
//-   Initial release
//- V1.1   28 Jan 2010
//-   If opposition point object already exists, update its value.

const DegreesToRadians = 0.017453292519943295769
const RadiansToDegrees = 57.295779513082320877
const HoursToRadians = 0.2617993878
const RadiansToHours = 3.8197186342

Function Normalise(x as double, vfrom as double, vto as double) As double
  dim norm,range as double
  
  norm = x
  range = vto-vfrom
  while norm<vfrom
    norm=norm+range
  wend
  while norm>vto
    norm=norm-range
  wend
  return norm
End Function

sub EclipticToEquatorial(lambda as double, beta as double, eps as double, byref ra as double, byref dec as double)
  dec=asin(sin(beta*DegreesToRadians)*cos(eps*DegreesToRadians)+cos(beta*DegreesToRadians)*sin(eps*DegreesToRadians)*sin(lambda*DegreesToRadians))*RadiansToDegrees
  ra = Normalise(atan2(sin(lambda*DegreesToRadians)*cos(eps*DegreesToRadians)-tan(beta*DegreesToRadians)*sin(eps*DegreesToRadians),cos(lambda*DegreesToRadians))*RadiansToHours,0.0,24.0)
end sub


sub EquatorialToEcliptic(ra as double,dec as double, eps as double, byref lambda as double, byref beta as double)
  lambda = atan2(sin(ra*HoursToRadians)*cos(eps*DegreesToRadians)+tan(dec*DegreesToRadians)*sin(eps*DegreesToRadians), cos(ra*HoursToRadians))*RadiansToDegrees
  beta=asin(sin(dec*DegreesToRadians)*cos(eps*DegreesToRadians)-cos(dec*DegreesToRadians)*sin(eps*DegreesToRadians)*sin(ra*HoursToRadians))*RadiansToDegrees
end sub

sub main()
  dim ra,dec,lambda,beta as double
  dim ob as APPlanObject, i as integer
  
  const eps = 23.4
  
  // Get current position of the sun
  SolarSystemPosition(ss_Sun,ra,dec)

  // Convert to ecliptic coordinates  
  EquatorialToEcliptic(ra,dec,eps,lambda,beta)
  
  // Opposition point is 180 deg along the ecliptic from the sun
  lambda=lambda+180.0
  beta = 0.0
  
  // Convert back to RA/Dec
  EclipticToEquatorial(lambda,beta,eps,ra,dec)
  
  // Check if the opposition point object is already there
  ob=nil
  for i = 1 to nObj
    if Obj(i).ID="OP" then
      ob=Obj(i)
      exit
    end if
  next
  
  if ob=nil then ob=NewObject
   
  ob.ID="OP"
  ob.Name="Opposition point"
  ob.RA=ra
  ob.Dec=dec
end sub

main