import jdcal
import numpy
import pyfits
import os
import datetime
import time
from array import array




u19700101=62135683200.0

ct = time.gmtime()
ct2 = time.strftime("%Y-%m-%d",ct)
ct2_1 = time.strftime("%Y.%m.%d",ct)
ct2_H = float(time.strftime("%H",ct))
ct2_M = float(time.strftime("%M",ct))
ct2_S = float(time.strftime("%S",ct))
ct3=time.time()
dt=datetime.datetime.strptime(ct2_1,'%Y.%m.%d')
julian=sum(jdcal.gcal2jd(dt.year,dt.month,dt.day))+(ct2_H+ct2_M/60.0+ct2_S/3600.0)/24.0



int_number_nx = 100
int_number_ny = 100
int_number_scan = int_number_nx*int_number_ny
int_size_spectrum = 4096

int_number = numpy.zeros((int_number_scan,),dtype=numpy.int)
int_scan = numpy.zeros((int_number_scan,),dtype=numpy.int)
str_datatype = numpy.zeros((int_number_scan,),dtype=numpy.dtype('S10'))
int_quality = numpy.zeros((int_number_scan,),dtype=numpy.int)
float_utobs = numpy.zeros((int_number_scan,),dtype=numpy.float64)
float_lstobs = numpy.zeros((int_number_scan,),dtype=numpy.float64)
float_azimuth = numpy.zeros((int_number_scan,),dtype=numpy.float64)
float_elevation = numpy.zeros((int_number_scan,),dtype=numpy.float64)
float_tsys = numpy.zeros((int_number_scan,),dtype=numpy.float64)
float_time = numpy.zeros((int_number_scan,),dtype=numpy.float64)
str_source = numpy.zeros((int_number_scan,),dtype=numpy.dtype('S20'))
float_lambda = numpy.zeros((int_number_scan,),dtype=numpy.float64)
float_latitude = numpy.zeros((int_number_scan,),dtype=numpy.float64)
float_off_lambda = numpy.zeros((int_number_scan,),dtype=numpy.float64)
float_off_latitude = numpy.zeros((int_number_scan,),dtype=numpy.float64)
float_epoch = numpy.zeros((int_number_scan,),dtype=numpy.float64)
str_line = numpy.zeros((int_number_scan,),dtype=numpy.dtype('S20'))
int_channels = numpy.zeros((int_number_scan,),dtype=numpy.int)
float_reference = numpy.zeros((int_number_scan,),dtype=numpy.float64)
float_freq_step = numpy.zeros((int_number_scan,),dtype=numpy.float64)
float_velo_step = numpy.zeros((int_number_scan,),dtype=numpy.float64)
float_velocity = numpy.zeros((int_number_scan,),dtype=numpy.float64)
float_frequency = numpy.zeros((int_number_scan,),dtype=numpy.float64)
float_beam_eff = numpy.zeros((int_number_scan,),dtype=numpy.float64)
float_pressure = numpy.zeros((int_number_scan,),dtype=numpy.float64)
float_ambient_t = numpy.zeros((int_number_scan,),dtype=numpy.float64)
float_spectrum = numpy.random.random((int_number_scan,int_size_spectrum,2,1,1))


for i in xrange(int_number_nx):
   for j in xrange(int_number_nx):
       ct = time.gmtime()
       ct2 = time.strftime("%Y-%m-%d",ct)
       ct2_1 = time.strftime("%Y.%m.%d",ct)
       ct2_H = float(time.strftime("%H",ct))
       ct2_M = float(time.strftime("%M",ct))
       ct2_S = float(time.strftime("%S",ct))
       ct3=time.time()
       dt=datetime.datetime.strptime(ct2_1,'%Y.%m.%d')
       julian=sum(jdcal.gcal2jd(dt.year,dt.month,dt.day))+(ct2_H+ct2_M/60.0+ct2_S/3600.0)/24.0
       int_number[i*int_number_nx+j] = 1
       int_scan[i*int_number_nx+j] = 1
       str_datatype[i*int_number_nx+j] = 'Line'
#        str_date[i*int_number_nx+j] = ct2
       int_quality[i*int_number_nx+j] = 1
       float_utobs[i*int_number_nx+j]= 1.0
       float_lstobs[i*int_number_nx+j]= 1.0
       float_azimuth[i*int_number_nx+j]= 1.0
       float_elevation[i*int_number_nx+j]= 1.0
       float_tsys[i*int_number_nx+j]= 1.0
       float_time[i*int_number_nx+j]= 1.0
       str_source[i*int_number_nx+j] = 'Orion'
       float_lambda[i*int_number_nx+j]=i*0.5/60.0
       float_latitude[i*int_number_nx+j]=j*0.5/60.0
       float_off_lambda[i*int_number_nx+j]=0.0
       float_off_latitude[i*int_number_nx+j]=0.0
       float_epoch[i*int_number_nx+j]=2000.0
       str_line[i*int_number_nx+j] = 'OH1665'
       int_channels[i*int_number_nx+j] = int_size_spectrum
       float_reference[i*int_number_nx+j] = 1.0
       float_freq_step[i*int_number_nx+j] =  1.0
       float_velo_step[i*int_number_nx+j] =  1.0
       float_velocity[i*int_number_nx+j] =  1.0
       float_frequency[i*int_number_nx+j] =  1.0
       float_beam_eff[i*int_number_nx+j] =  1.0
       float_pressure[i*int_number_nx+j]=1.0e5
       float_ambient_t[i*int_number_nx+j]=1.0
       #time.sleep(5.1)

os.system('rm -rf gen_scantable.fits')
hdulist = pyfits.open('test_scantable.fits')
hdu1 = hdulist[1]
tempphdu = hdulist[0]
phheader=pyfits.Header()
phdu = pyfits.PrimaryHDU(header=phheader)

column1_data = pyfits.Column(name='NUMBER',format='1I',array=int_number)
column2_data = pyfits.Column(name='SCAN',format='1I',array=int_scan)
column3_data = pyfits.Column(name='DATATYPE',format='10A',array=str_datatype)
column4_data = pyfits.Column(name='QUALITY',format='1I',array=int_quality)
column5_data = pyfits.Column(name='UTOBS',format='1E',unit='s',array=float_utobs)
column6_data = pyfits.Column(name='LSTOBS',format='1E',unit='s',array=float_lstobs)
column7_data = pyfits.Column(name='AZIMUTH',format='1E',unit='degree',array=float_azimuth)
column8_data = pyfits.Column(name='ELEVATION',format='1E',unit='degree',array=float_elevation)
column9_data = pyfits.Column(name='TSYS',format='1E',unit='K',array=float_tsys)
column10_data = pyfits.Column(name='TIME',format='1E',unit='s',array=float_time)
column11_data = pyfits.Column(name='SOURCE',format='20A',array=str_source)
column12_data = pyfits.Column(name='LAMBDA',format='1E',unit='degree',array=float_lambda)
column13_data = pyfits.Column(name='LATITUDE',format='1E',unit='degree',array=float_latitude)
column14_data = pyfits.Column(name='OFF_LAMBDA',format='1E',unit='degree',array=float_off_lambda)
column15_data = pyfits.Column(name='OFF_LATITUDE',format='1E',unit='degree',array=float_off_latitude)
column16_data = pyfits.Column(name='EPOCH',format='1E',unit='year',array=float_epoch)
column17_data = pyfits.Column(name='LINE',format='10A',array=str_line)
column18_data = pyfits.Column(name='CHANNELS',format='1I',array=int_channels)
column19_data = pyfits.Column(name='FREQ_STEP',format='1E',unit='MHz',array=float_freq_step)
column20_data = pyfits.Column(name='VELO_STEP',format='1E',unit='km/s',array=float_velo_step)
column21_data = pyfits.Column(name='VELOCITY',format='1E',unit='km/s',array=float_velocity)
column22_data = pyfits.Column(name='FREQUENCY',format='1E',unit='MHz',array=float_frequency)
column23_data = pyfits.Column(name='BEAM_EFF',format='1E',array=float_beam_eff)
column24_data = pyfits.Column(name='PRESSURE',format='1D',unit='Pa',array=float_pressure)
column25_data = pyfits.Column(name='AMBIENT_T',format='1D',unit='K',array=float_ambient_t)


column26_data = pyfits.Column(name='DATA',format='8192E',dim='(4096,2,1,1)',unit='K',array=float_spectrum)



table_hdu = pyfits.new_table([column1_data,column2_data,column3_data,column4_data,column5_data,column6_data,column7_data,column8_data,column9_data,column10_data,column11_data,column12_data,column13_data,column14_data,column15_data,column16_data,column17_data,column18_data,column19_data,column20_data,column21_data,column22_data,column23_data,column24_data,column25_data,column26_data])


table_hdu.header.insert('TTYPE1',('EXTNAME','SINGLE DISH','name of this binary table extension'))
table_hdu.header.insert('TTYPE1',('NMATRIX',1,'Number of DATA arrays'))
table_hdu.header.insert('TTYPE1',('OBSERVER','QL','Observer name(s)'))
table_hdu.header.insert('TTYPE1',('PROJID','1-1-1','Project name'))
table_hdu.header.insert('TTYPE1',('TELESCOP','FAST','Telescope name'))
table_hdu.header.insert('TTYPE1',('OBSGEO-X',-1668557.9202,'[m] Antenna ITRF X-coordinate'))
table_hdu.header.insert('TTYPE1',('OBSGEO-Y',5506865.9640,'[m] Antenna ITRF X-coordinate'))
table_hdu.header.insert('TTYPE1',('OBSGEO-Z',2744946.0867,'[m] Antenna ITRF X-coordinate'))
table_hdu.header.insert('TTYPE1',('CTYPE1','FREQ','DATA array axis 1: frequency in Hz'))
table_hdu.header.insert('TTYPE1',('CTYPE2','STOKES','DATA array axis 2: polarization code'))
table_hdu.header.insert('TTYPE1',('CRPIX2',1.0,'Polarization code reference pixel'))
table_hdu.header.insert('TTYPE1',('CRVAL2',-5.0,'Polarization code at reference pixel (XX)'))
table_hdu.header.insert('TTYPE1',('CDELT2',-1.0,'Polarization code axis increment'))
table_hdu.header.insert('TTYPE1',('CTYPE3','RA','DATA array axis 3 (degenerate):RA (mid-int)'))
table_hdu.header.insert('TTYPE1',('CRPIX3',1.0,'RA reference pixel'))
table_hdu.header.insert('TTYPE1',('CDELT3',-1.0,'RA axis increment'))
table_hdu.header.insert('TTYPE1',('CTYPE4','DEC','DATA array axis 4 (degenerate):Dec (mid-int)'))
table_hdu.header.insert('TTYPE1',('CRPIX4',1.0,'Dec reference pixel'))
table_hdu.header.insert('TTYPE1',('CDELT4',1.0,'Dec axis increment'))
table_hdu.header.insert('TTYPE1',('SPECSYS','LSRK','Doppler reference frame (transformed)'))
table_hdu.header.insert('TTYPE1',('SSYSOBS','TOPOCENT','Doppler reference frame of observation'))
table_hdu.header.insert('TTYPE1',('EQUINOX',2000.0,'Equinox of equatorial coordinates'))
table_hdu.header.insert('TTYPE1',('RADESYS','FK5','Equatorial coordinate frame'))




hdulist2 = pyfits.HDUList([phdu,table_hdu])
hdulist2.writeto('gen_scantable.fits')

gen_scantable.zip