import numpy as np, sys
#
print '\n time-estimator-IRAM30m-nika-gismo v 2014.03.28'
print '\n Based on the Billot et al. (2014) time estimator guide. \n'
#
if len(sys.argv) != 17:
    print ' USAGE:\n'
    print '   python time-estimator-IRAM30m-nika-gismo.py --camera "GISMO" --rms 0.5 --pwv 2 --elevation 30 --Xsize 4 --Ysize 6 --filtering 1.2 --overheads 2.0'
    print '\n'
    print ' OPTIONS:\n'
    print '        camera =>  IRAM 30m continuum camera. Options are GISMO, NIKA1mm, NIKA2mm.'
    print '           rms =>  Any value above the confussion limit.'
    print '           pwv =>  Precipitable water vapor. Possible values 2, 4, and 7 mm.'
    print '     elevation =>  Values from 20 to 83 degrees.'
    print '   Xsize Ysize =>  Mapping area. Xsize and Ysize must be equal (point sources) or larger than the Field of View.'
    print '     filtering =>  Data filtering scheme. Values from 1 to 4.'
    print '     overheads =>  Telescope/observing overheads (calibration, pointing, focus). Values from 1.6 to 2.6.\n'
    sys.exit()
#
#
for i in range(8):
#
   if sys.argv[i*2+1] == ("--camera"):
      camera = sys.argv[i*2+2]
#
   if sys.argv[i*2+1] == ("--rms"):
      rms = float(sys.argv[i*2+2])
#
   if sys.argv[i*2+1] == ("--pwv"):
      if float(sys.argv[i*2+2]) < 3:
         pwv=2
      if 3 <= float(sys.argv[i*2+2]) < 5.5:
          pwv=4
      if float(sys.argv[i*2+2]) >= 5.5:
          pwv=7
#
   if sys.argv[i*2+1] == ("--elevation"):
      elevation = float(sys.argv[i*2+2])
#
   if sys.argv[i*2+1] == ("--Xsize"):
      Xsize = float(sys.argv[i*2+2])
#
   if sys.argv[i*2+1] == ("--Ysize"):
      Ysize = float(sys.argv[i*2+2])
#
   if sys.argv[i*2+1] == ("--overheads"):
      overheads = float(sys.argv[i*2+2])
#
   if sys.argv[i*2+1] == ("--filtering"):
      filtering = float(sys.argv[i*2+2])
#
#
#
if camera == 'GISMO':
#
    Np=100.                     #Number of pixels
    Sp=13.75                    #Pixel size (in arcsec)
    HPBW=16.7                   #Half-power beamwidth (in arcsec)
    Wavelength=1200.            #Central wavelength in microns
    NEFD=14.                    #Noise ==uivalent flux density
#
    if pwv == 2:
       tau=0.06
    if pwv == 4:
       tau=0.11
    if pwv == 7:
       tau=0.19  
#
    if Xsize <= 1.:
        Xsize=1.8
        Xpoint=1
    else:
        Xpoint=0
#
    if Ysize <= 1.8:
        Ysize=3.7
        Ypoint=1
    else:
        Ypoint=0
#
#
#
if camera == 'NIKA1mm':
#
    Np=136.                     #Number of pixels
    Sp=6.8                      #Pixel size (in arcsec)
    HPBW=12.                    #Half-power beamwidth (in arcsec)
    Wavelength=1250.            #Central wavelength in microns
    NEFD=35.                    #Noise equivalent flux density
#
    if pwv == 2:
        tau=0.22
    if pwv == 4:
        tau=0.36
    if pwv == 7:
        tau=0.57  
#
    if Xsize <= 1.:
        Xsize=1.8
        Xpoint=1
    else:
        Xpoint=0
#
    if Ysize <= 1.:
        Ysize=1.8
        Ypoint=1
    else:
        Ypoint=0
#
#
#
if camera == 'NIKA2mm':
#
    Np=114.                     #Number of pixels
    Sp=9.6                      #Pixel size (in arcsec)
    HPBW=17.5                   #Half-power beamwidth (in arcsec)
    Wavelength=2000.            #Central wavelength in microns
    NEFD=14.                    #Noise equivalent flux density
#
    if pwv == 2:
        tau=0.06
    if pwv == 4:
        tau=0.11
    if pwv == 7:
        tau=0.19  
#
    if Xsize <= 1.0:
        Xsize=2.0
        Xpoint=1
    else:
        Xpoint=0
#
    if Ysize <= 1.0:
        Ysize=2.0
        Ypoint=1
    else:
        Ypoint=0
#
#
transmission = np.e**(-tau/np.sin(elevation*np.pi/180.))
#
if Xpoint+Ypoint != 2:
   time=(NEFD*filtering/transmission/rms)**2*(1.+Xsize*60.*Ysize*60./Np/Sp**2)*overheads/3600.
   print '        camera => ', camera
   print '           rms => ', rms, ' mJy/beam'
   print '           pwv => ', pwv, ' mm'
   print '     elevation => ', elevation, ' deg'
   print '   Xsize Ysize => ', Xsize, Ysize, ' arcmin'
   print '     overheads => ', overheads
   print '     filtering => ', filtering
   print '                 '
else:
   time=(NEFD*filtering/transmission/rms)**2*overheads/3600.
   print '        camera => ', camera
   print '           rms => ', rms, ' mJy/beam'
   print '           pwv => ', pwv, ' mm'
   print '     elevation => ', elevation, ' deg'
   print '   Xsize Ysize => ', 'Point source'
   print '     overheads => ', overheads
   print '     filtering => ', filtering
   print '                 '
   #
if time >= 1.:
   time_str=str(round(time,1))+' hours'
#
if time <= 1. and time >= 1./60.:
   time_str=str(round(time*60,1))+' minutes'
#
if time <= 1./60. and time >= 0.1/3600.: 
   time_str=str(round(time*3600,1))+' seconds'
#
if time <= 0.1/3600.:
   time_str='0.1 seconds'
#
print '          time => '+time_str+' \n'
#
print 'Please include the following text into your proposal:\n'
#
print 'The total observing time using '+camera+' to map a region of '+str(Xsize)+'x'+str(Ysize)+' arcminutes to reach an rms of '+str(rms)+' mJy/beam, assuming '+str(pwv)+' mm pwv, '+str(elevation)+' deg elevation, Ffilter='+str(filtering)+', Foverhead='+str(overheads)+', was estimated to be *'+time_str+'*, using the time estimator v 2014.03.27. \n '