text: #! /usr/bin/env python
  
  import time
  import numpy
  import ephem
  
  
  def sunradec(utime):
      """Return the [RA,Dec] (in radians) of the Moon for a given unix time"""
      sun = ephem.Sun()
      utime = time.localtime(utime) # get the time tuple
      sun.compute((utime[0],utime[1],utime[2],utime[3],utime[4],utime[5]))
      return [sun.ra/1.0,sun.dec/1.0]
  
  # Vector operations
  
  def cross(v2,v1):
      """Calculate the cross product of 2 3D vectors"""
      vec = numpy.array([0,0,0],numpy.float)
      vec[0] = v2[1]*v1[2]-v2[2]*v1[1]
      vec[1] = v2[2]*v1[0]-v2[0]*v1[2]
      vec[2] = v2[0]*v1[1]-v2[1]*v1[0]
      return vec
  
  def vecnorm(vec):
      """Normalise a vector."""
      mag = numpy.sqrt(sum(vec**2))
      return vec/mag
  
  def radec2vec(ra,dec):
      """Convert RA/Dec angle (in radians) to a vector"""
  
      v1 = numpy.cos(dec)*numpy.cos(ra)
      v2 = numpy.cos(dec)*numpy.sin(ra)
      v3 = numpy.sin(dec)
  
      return numpy.array([v1,v2,v3])
  
  def ICSdateconv(date):
      """Convert the date format used in the ICS to standard UNIX time"""
      x = date.replace("/"," ").replace("-"," ").replace(":"," ").split()
      base = time.mktime((int(x[0]),1,0,0,0,0,0,0,0))
      return base + ((int(x[1]))*86400 + (float(x[2])) * 3600 + float(x[3])*60 + float (x[4]))
  
  
  def optimum_roll(ra,dec,utime):
      """Calculate the optimum Swift Roll angle (in radians) for a given Ra, Dec and Unix Time"""
      # Calcuate the vector for the Sun position
      sun = sunradec(utime)
      vSun = radec2vec(sun[0],sun[1])
  
      # define the target vector
      vT = radec2vec(ra,dec)
  
      # get cross product of Sun vector and target vector (vSun x vT) :  result is vector Y
  
      vY = cross(vT,vSun)
      vnY = vecnorm(vY) # Normalise the Y vector
  
      # get cross product of normalized Y and target vector: result is vector Z
  
      vZ = cross(vnY,vT)
      vnZ = vecnorm(vZ)
  
      # newroll = -atan(vnY[2]/vnZ[2]), to get proper newroll use function atan2
      # JAK Note: For some reason this doesn't work - so I got rid of the sign and it did.
  
      if (vnY[2] != 0 and vnZ[2] != 0):
          newroll = numpy.arctan2(vnY[2],vnZ[2])
          newroll = newroll/dtor
      else:
          # roll is not uniquely defined, arbitrarily pick 0.0 or 180.0
          newroll = 0
          if (vSun[0]*(-cos(newroll)*sin(dec)*cos(ra) - sin(newroll)*sin(ra))
              + vSun[1]*(-cos(newroll)*sin(dec)*sin(ra) + sin(newroll)*cos(ra))
              + vSun[2]*cos(newroll)*cos(dec) < 0.0):
              # 0.0 would put the sun vector in -Z, so use 180.0 deg
              newroll = 180.0
  
      if (newroll < 0.0):
          newroll += 360.0
      return newroll*dtor
  
  
  dtor = numpy.pi/180.0
  
  lt = time.gmtime()
  year = lt[0]
  day = lt[7]
  
  #ra = 266*dtor
  #dec = -29*dtor
  
  rollcons = 9.5
  suncons = 45
  
  def table(ra,dec,year=year,startday=0,endday=265):
      """ Produce a table of approximate roll ranges for Swift observation of a target 
          with given position.
          
          Parameters:
          
          ra, dec : float
             position in decimal degrees (J2000)
          year : int (optional)
             year for the table (default = current year)     
      """
      print("Roll angle ranges for swift observation of target at\n RA,Dec=(%.2f,%.2f) year=%i\n"%
         (ra,dec,year))
      ra = ra*dtor
      dec = dec*dtor
      print("%3s %7s %7s %7s %7s %7s %7s"%("Day","SunRA","SunDec","SunDist","Roll","MinRoll","MaxRoll"))
      for i in range(startday,endday+1):
          now = ICSdateconv("%s-%03d-00:00:00"%(year,i))
          sra,sdec = sunradec(now)
          oroll = optimum_roll(ra,dec,now)
          sundist = ephem.separation([ra,dec],[sra,sdec])
          roll_range = rollcons/numpy.sin(sundist);
          if sundist/dtor > suncons:
              #print(f"{i:3d} {sra/dtor:7.3f} {sdec/dtor:7.3f} {sundist/dtor:7.3f} 
              #{oroll/dtor:7.3f} {oroll/dtor-roll_range:7.3f} {oroll/dtor+roll_range:7.3f}")
              print("%3d %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f"%(i,sra/dtor,sdec/dtor,
                sundist/dtor,oroll/dtor,oroll/dtor-roll_range,oroll/dtor+roll_range))
  
  
  
  