'''PlateSystem
--- 7 June 2006 KIDS with eye on second  decimal in field
----            plots produce contours of fieldsize_ra and fieldsize_dec
--- 1 June 2006 KIDS adaptations
--- for KIDS different overlap in RA and DEC   1 june 2006
--- SMOOTHING DISABLED- not needed FOR SUCH SMALL OVERLAP
--- overlap defined at BOTTOM of fields -> more fields/ more overlap
Define a platesystem, a collection fields ("plates") of a constant size
covering the entire sky.
Edwin Valentijn -
---- 7 March 2004
---- version with smooth grid - more advanced than PAE
'''
__version__ = '@(#)$Revision: 1.3 $'

import math

def sum(args):
   '''The sum function is ONLY built in in Python 2.3 and newer'''
   import operator
   return reduce(operator.add, args, 0)

class PlateSystem:
    '''Define a plate system

    The plate system is defined by a set of fields covering the sky.
    Each field has a size in ra and dec, a centre and a requested amount of overlap with
    neighbouring fields in both ra and dec and a minimun allowed overlap in ra
    The fields are arranged in strips perpendicular to the meredian.
    This means that with increasing declination there are less fields per strip.
    Script optimizes for continuity form one dec strip to the next;
    therfor the positions in various dec stips are dependent on each other.
    '''
    def __init__(self, fieldsizera=1.0127, fieldsizedec=1.02451, \
                 overlapra=1.186,overlapdec=4, efflimit=0.96):
        '''fieldsize - the size of the output fields in degrees
        overlap      -  % fraction of overlap between fields
        efflimit     -  % fraction minimum overlap at smallest dec in the field
         '''
        self.fieldsizera    = fieldsizera                        # 1.0
        self.fieldsizedec   = fieldsizedec                       # 1.0
        self.efflimit       = efflimit                           # 0.0
        self.overlapra      = 1. +overlapra/100.                 # 1% -> 1.01
        self.overlapdec     = 1. +overlapdec/100.                # 5% -> 1.05
        self.distancera     = fieldsizera /  self.overlapra      # 1/1.01 ->0.99
        self.distancedec    = fieldsizedec / self.overlapdec     # 1/1.05 ->0.95
        print "requested overlap ra=", overlapra,"%  requested overlap dec=", overlapdec, "%"
# make dec centers grid first
        numberstrips        = int(90.0 / self.distancedec)
        dec                 = range (numberstrips)
        for j in range (0, numberstrips):
#            dec[j]    = j * 90.0/numberstrips
# 1 June 2006 KIDS start at dec=0.5 instead of 0.0
             dec[j]    = 0.5+(j * 90.0/numberstrips)
# 1 June 2006 KIDS remove polar field (now dec start at 0.5)            
#        self.numberstrips   = numberstrips +1                       # add polar field
#        dec.append(89.9)                                            # add polar
        self.numberstrips   = numberstrips
        self.dec            = dec                                   # key grid of stips in de
#        print self.dec
        d = dec[1]-dec[0]
        print "fieldsize dec =", fieldsizedec,"degrees,  fieldsizera=", \
              fieldsizera, "degrees"
        print "achieved distance fieldcenters in dec:", d, "degrees"
        print "achieved overlap in dec:",  100.*(fieldsizedec-d)/d, "%"
        self.n,self.e       = self.get_numberra_centers(0)
#
#  solution for numberoffields in RA for each dec strip
#
    def get_numberra_centers(self,log):    # determine number of fields for each dec
        '''
        log   - print log to screen '''
        nold     = 400
        dold     = 0
        dec      = self.dec 
        n =  range (self.numberstrips)
        e =  range (self.numberstrips)
        for i in range (0, self.numberstrips):
#         smallcircle      = 360.0 * math.cos( self.dec[i]* math.pi / 180.)
# 1 June 2006 improve take smallcircle at bottom of field !!!
           smallcircle      = 360.0 * math.cos( (self.dec[i]-self.fieldsizedec/2.0)* math.pi / 180.)
           numberperstrip   = int(round(smallcircle / self.distancera))    
           effectiveoverlap = 100*( (self.fieldsizera/(smallcircle/numberperstrip)) -1.)
           if log:
             print "--------------------------------------------------------------------"
             print "rawoverlap-ra in %=", "%5.2f" %effectiveoverlap
                                            # finetune on grid 1
           while effectiveoverlap < self.efflimit :       # handles both negative and limit on overlap
                numberperstrip  +=1
                effectiveoverlap = 100*( (self.fieldsizera/(smallcircle/numberperstrip)) -1.)
                if log: print "plusone overlap"
                                            # finetune on grid 2
                                            # smooth out the jumps in strips
                                            # CREATES SMOOTH GRID  but  not for pole
# 1 June 2006 smooth removed for small overlap
          # if (nold-numberperstrip)<dold and dec[i]<84 :
          #      numberperstrip   = nold - dold
          #      effectiveoverlap = 100*( (self.fieldsize/(smallcircle/numberperstrip)) -1.)
          #      if log:
          #         print "plusone smooth"
           if log:
                epix=16000.*effectiveoverlap/100.
                print "                                                         delta",(nold -numberperstrip)
                print i+1,\
                      "dec=", "%4.2f" %self.dec[i], \
                      "small circle=" ,"%5.2f" %smallcircle,\
                      "numberperstrip=",numberperstrip, \
                      "overlap-ra%=", "%5.2f" %effectiveoverlap,\
                      "%3.0f" %epix,"px"
           if i > 0 :
                dold = nold-numberperstrip 
           nold = numberperstrip
           n[i] = numberperstrip
           e[i] = effectiveoverlap
        if log:
           print "PlateSystem: total number found=",sum(n)
        self.n = n
        self.e = e
        return n,e
#
    def nearest_dec(self,decin):
        ''' get nearest fieldcenter dec  for any dec'''
        dec  = self.dec
        decin = abs (decin)                 # allow pass negative dec
        diff = range( len( dec))
        for i in range (0,len( dec)):
            diff[i] = abs( decin-dec[i])
        j = diff.index(min(diff))
        return dec[j],j
#        
    def get_ra(self,decin):                 # get seq of ra's for given decin
        ''' get seq of ra's for a single decin'''
        decplate,j = self.nearest_dec(decin)# get numberperstrip and index
        n = self.n
        ras             = range (n[j])        
        offsetra        = 360./n[j]
        for i in range (0,n[j]):
            ras[i] = i * offsetra
        return ras,decplate
#    
    def get_field_centre(self,rain,decin):        
        ''' get nearest ra fieldcenters, needs dec too'''
        if decin < 0:
            decsign = -1
        else:
            decsign = 1
        dec = abs(decin)
        if dec > 90.0:
            print" dec >90 we can't go there"
            dec = 90.
        ras, decplate = self.get_ra(dec)
        diff = range(len(ras))       
        for i in range(0, len(ras)):
            diff[i] = abs(rain-ras[i])
        j = diff.index(min(diff))
        rasplate = ras[j]
        return rasplate, decsign*decplate
#    
    def print_grid(self):
        ''' print grid properties in file'''
        print "Version for parents of KIDS " 
        print "sequence     Dec      number       overlap-ra"
        print "         [degrees]                 [min%] min[px] +top[px]"       
        for i in range (0, self.numberstrips):
           top = math.cos((self.dec[i] - self.fieldsizedec/2.0) * math.pi / 180.)-\
                 math.cos((self.dec[i] + self.fieldsizedec/2.0) * math.pi / 180.)
           top*=16000.
           print "%8i" %(i+1), "%9.2f" %self.dec[i],"%9i" %self.n[i], \
                 "%9.1f" %self.e[i], "%5.0f " %(160*self.e[i]),"%5.0f "%top
        print  "PlateSystem: total number of fields=", sum(self.n)
        return
# plot
# input parameters - range on sky grid -for plot only
    def plot_it(self,rastart=-1.,raend=360.,decstart=0.,decend=90.,title=' '):
        ''' plots the grid in range
        can be used as survey preparation tool
        selects fields within rectangle on sky
        returns ra,dec of selected fields in self
        autoscaling implemented                   ''' 

        import dislin

        rastart = float(rastart)
        decstart= float(decstart)
        raend   = float(raend)
        decend  = float(decend)
        selected_fields_ra = range(50000)
        selected_fields_dec= range(50000)
        if abs(decend-decstart)<1.0 :
            decend  +=0.5
            decstart-=0.5
        if abs(raend-rastart)<1.0:
            rastart -=0.5
            raend   +=0.5
        symsize =   12
        symtype =    8
        xlen    = 2200
#
# outfile (metafile)format ; befor disini 'post'=ps 'pdf'
        #dislin.metafl('png')      # outfile format ; befor disini
                                  # gif, pdf eps, png, tiff 
        dislin.metafl('cons')        # on screen  full screen size
        #dislin.metafl('xwin')        #    screen small window
        #dislin.setfil('C:\myPy\skygridpole.tiff')   # default is dislin
        dislin.disini ()
        dislin.hwfont ()             # standard hardware fonts
        #?dislin.page(6400,200)         # ?? does it do anything ??
        dislin.pagera ()             # plot border
        dislin.axspos (450, 1800)    # lower left corner
        dislin.axslen (xlen, 1400)   # length height of plot
        dislin.name   ('RA degrees' ,'X')
        dislin.name   ('Dec degrees','Y')
        dislin.digits (-1, 'XY')     # -1 integer
        dislin.ticks  (3, 'X')       # number of ticks between labels
        dislin.ticks  (10, 'Y')      # number of ticks between labels
        dislin.texmod ('ON')
        dislin.titlin ('OmegaCam \it Sky Grid  1.0 degree', 1)
        ticksx = 20
        if abs(raend-rastart)<100: ticksx = 10     # autoscaling very premature
        if abs(raend-rastart)<10: ticksx = 1      
        ticksy = 10.
        if abs(decend-decstart) < 15: ticksy = 1.
        dislin.graf   (rastart,raend, 10*int(rastart/10), ticksx,\
#                       decstart, decend, 10*int((decstart-10.0)/10), ticksy) # axi
                        decstart,decend,decstart,ticksy)
        dislin.title  ()             # print title defined in titlin - up to 4 lines
        dislin.color  ('blue')
        dislin.hsymbl (symsize)           # symbol size
#
        y    = self.dec
# 1 June 2006 start at dec =0.5 instead 0.0        
#        ymin = range (len(y)-1)         # append negative dec s
#        for ii in range(len(ymin)): ymin[ii] = -1 * y[ii+1]
        ymin = range (len(y))          # append negative dec s
        for ii in range(len(ymin)): ymin[ii] = -1 * y[ii]
        ymin.extend(y)
        ymin.sort()
        y  =  ymin                        # append done
        fr  = self.fieldsizera
        fd  = self.fieldsizedec
        count = 0 
        dc = 1400/abs(decend-decstart) #pix/degree for plot coord (instead of user)
        rc = xlen/abs(raend-rastart)   #pix/degree
        for i in range (0, len(y)):
            y1 = y[i]
            cosdec= math.cos( y1 * math.pi / 180.)
            x,de   = self.get_ra(y1)
            for j in range (0, len(x)):
                if x[j]> rastart and x[j] < raend and y1 >decstart and y1<decend:
                    count += 1 
                    selected_fields_ra[count-1] = x[j]
                    selected_fields_dec[count-1]= y1
                    dislin.rlsymb (symtype,x[j],y1)   #  symbols in user coord
                    dislin.color ('green')
                    dislin.rectan (int(450+( (x[j]-(rastart+(0.5*fr/cosdec))) *rc)),\
                                   int(1800-((0.5*fd +y1-decstart)*dc)), \
                                   int(fr*rc/cosdec),\
                                   int(fd*dc))  # rectangles in plotcoord
                    dislin.color  ('blue')
        self.selected_field_ra  = selected_fields_ra[0:count]
        self.selected_field_dec = selected_fields_dec[0:count]
        e = self.e
        e = e[0:len(e)-1]             # skip polar field for overlap computation
        x = 'PlateSystem: # pointings=  '+str(sum(self.n))+ \
            ' average overlap ra=  '+str(sum(e)/len(e))[0:4]+'%'\
           ' selected pointings='+str(count)+' '+title
        print x
        dislin.titlin(x,3)
        dislin.titlin('KIDS optimized minimum overlap, different in RA and Dec',2)
        dislin.titlin(title,4)
        dislin.title  ()
        dislin.dash   ()
        dislin.xaxgit ()
        dislin.disfin ()
        return 
#    
# ways of invoking the class
#                              
if __name__ == '__main__':
    grid = PlateSystem()
#    a = grid.get_numberra_centers(1)      # makes full solution - print log
    a = grid.get_numberra_centers(0)       #  full solution without log
#   nearest dec for free input dec
#    for i in range(1,9):
#        print "dec fieldccenter", 9.*i, grid.nearest_dec(9.*i)
#        print "ra- dec nearest ra =", 9.*i,'dec=',9.*i ,grid.get_field_centre(9.*i,9.*i)
#
# print the key table of the plate grid
    a= grid.print_grid()    
# plotit        
    b= grid.plot_it()

# survey preparation - list of coordinates: example of preparing KIDS survey"
    b = grid.plot_it(150,238,-5,3,'KIDS -NORTH1')  # N 
 #   b = grid.plot_it(225,238,-5,-3,'KIDS -NORTH1-not')  # N 1 not
#    for i in range(0,len(grid.selected_field_ra)):
#	       print i+1,grid.selected_field_ra[i],grid.selected_field_dec[i]
#
    b = grid.plot_it(130,137,-3,-1,'KIDS -NORTH2') # N 2
#    for i in range(0,len(grid.selected_field_ra)):
#	       print i+1,grid.selected_field_ra[i],grid.selected_field_dec[i]
#

#    b = grid.plot_it(330,360,-38,-25,'KIDS -SOUTH+SPT-part1')   # -35 without SPT
#    b = grid.plot_it(-0.5,52.5,-38,-25,'KIDS -SOUTH+SPT-part2')
    b = grid.plot_it(239.50,242.5,15,19,'VESUVIO-Hercules')
#    b = grid.plot_it(41,56.25,-58,-44,'VESUVIO-Horologium')
#    b = grid.plot_it(0,360,89,90,'Pole')
#    b = grid.plot_it(28,30,28,30,'example')
