#!/usr/bin/env python

import sys
import math
import numpy as np
#import pdb

"""
Return lits with possible "ncomb" that can be used together to define reduced
grid, according to input resolution.

13 Jun 2013 - Ph. Le Sager - v0

"""

# ------- Err
def wrong_input(mess):
    print mess
    print("Usage: %s RESOLUTION_LONG RESOLUTION_LAT" % sys.argv[0])
    return

# -------
def replace_multi(mylist, reflist):
    """
    Return lists formed with input list + results of
    find_multiple applied to last element of input list
    (Recursive)
    """
    val=mylist[-1]
    a=reflist.index(val)
    zelists=[]
    b=[f for f in reflist[a+1:] if not f%val] # find_multiple
    for i,j in enumerate(b):

        # skip some (they would be subset of another)
        check=[k for k in b[:i] if j%k==0]
        if check: continue
        
        c=mylist[:]
        c.append(j)
        newlists=replace_multi(c,reflist)
        if newlists:
            for k in newlists: zelists.append(k)
        else:
            zelists.append(c)
    return zelists

# ------- Input
total = len(sys.argv)
if total < 3 :
    dummy=wrong_input("Two arguments required!")
    sys.exit(1)

resolution=float(sys.argv[1])

# May need different check for some non-integer resolution (1/3 for example):
if 360.%resolution:
    dummy=wrong_input("Argument must divide 360!")
    sys.exit(1)
    
nbox=int(360/resolution)
print ""
print "Resolution  = %s" % resolution
print "Nb of boxes = %s" % nbox


# ------- Possible sequences of multiples

divisors=[f for f in range(2,nbox+1) if not nbox%f] # List of divisors

alllists=[]

for i in range(len(divisors)):
    l=divisors[i:i+1]

    # limit to "prime number" in the list
    cl=[k for k in divisors[0:i] if not l[0]%k]
    if cl: continue

    newlists=replace_multi(l, divisors)

    for k in newlists: alllists.append(k)

if alllists:
    print """
    You can use any of the following %s sequences (or any subset of them) to
    define a ReducedGrid. But you cannot mix these sequences.
    """ % len(alllists)
    
    for k in alllists: print k


# ------- Choice and corresponding grid size

Re=6371. # Earth radius in km
pi=math.pi

lat_res=float(sys.argv[2])

choice=raw_input("""
Enter comma-separated list of 'ncomb', starting from pole.
(for example: 40,8,8,4,4,4,2). Leave empty to exit.
> """)

if choice == '' : sys.exit(0)
choice=choice.split(',')
choice=np.array(choice,dtype=int)
nbands=len(choice)

lats=90.-(np.arange(nbands)+0.5)*lat_res # at center of the grid boxes

ndx = 2*pi*Re*np.cos(lats*pi/180.)*choice/nbox # box width in km
odx = 2*pi*Re*np.cos(lats*pi/180.)/nbox # old box width in km

print ""
print " band #".rjust(8), "Lat.".rjust(7), "Old Width, New Width [km]"
for i,k in enumerate(ndx):
    print repr(i).rjust(8), ' {0:6.2f}'.format(lats[i]), \
          '{0:9.1f}'.format(odx[i]), '{0:10.1f}'.format(k)