import numpy
msfile='simple_obs.ms'
modelimage='glmtest.modim'
dir0 = me.direction('J2000','17h04m13.0','-42d19m58.0');
fbands         = ['Band3'];                # Names to give bands
freqs          = ['86.0GHz'];
nchan          = 8;                        # Number of spectral channels
dfreq0         = '2.0GHz';      # Channel width
refdates    = ['2012/05/21/05:05:13']; # Reference date/time
itime          = '10s'

###some fake calibrator at 
dircal=me.direction('J2000','17h30m13.0','-45d05m58.0');

cl.addcomponent(dir=dircal, flux=0.2, freq='86GHz')
cl.rename('mycal.cl')
cl.done()


##### define antennas
stnn = 50;
stnx = numpy.array([27789.81, 27806.26, 27801.31, 27818.16, 27848.62,
         27844.52, 27839.69, 27833.36, 27814.26, 27803.30,
         27861.26, 27860.31, 27824.72, 27840.40, 27834.85,
         27803.74, 27817.20, 27782.36, 27773.43, 27762.59,
         27780.13, 27780.58, 27764.71, 27749.51, 27739.30,
         27736.72, 27730.34, 27767.25, 27753.22, 27755.38,
         27728.02, 27720.42, 27717.98, 27746.80, 27797.54,
         27808.00, 27772.10, 27742.81, 27759.05, 27780.82,
         27800.04, 27825.62, 27840.30, 27743.11, 27759.60,
         27786.62, 27803.34, 27821.50, 27853.47, 27808.20]);

stny = numpy.array([53079.62, 53082.62, 53100.27, 53068.35, 53047.47,
        53070.45, 53091.26, 53104.90, 53128.29, 53138.64,
        53055.91, 53071.14, 53089.41, 53118.33, 53132.37,
        53158.27, 53107.16, 53066.47, 53081.30, 53069.82,
        53100.26, 53139.65, 53124.97, 53108.03, 53087.90,
        53060.58, 53046.95, 53147.06, 53136.80, 53094.18,
        53097.99, 53073.19, 53038.02, 53074.57, 53066.13,
        53045.89, 53045.88, 53027.71, 53021.54, 53014.65,
        53013.74, 53028.83, 53025.46, 53012.46, 53005.88,
        53030.28, 52997.74, 53003.23, 53017.97, 53027.93]);

stnz =numpy.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
       0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
       0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
       0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
       0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]);

stnd =numpy.array([12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0,
       12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0,
       12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0,
       12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0,
       12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0]);


sm.open(msfile)

posalma=me.observatory('ALMA')
n=stnn
ax=sum(stnx)/n;
x= stnx-ax; 
ay=sum(stny)/n;
y= stny-ay;
az=sum(stnz)/n;
z= stnz-az;
diam=stnd;

antnames=[]
for k in range(0,50) :
    antnames.append('ALMA%02d'%k)

###
### set the array configuration
sm.setconfig(telescopename='ALMA', x=x.tolist(), y=y.tolist(), z=z.tolist(), 
              dishdiameter=diam.tolist(), 
              mount=['alt-az'], antname=antnames,
              coordsystem='local', referencelocation=posalma);
###set spectralwindow
sm.setspwindow(spwname='band3', freq='86GHz', deltafreq='2GHz', 
                    freqresolution='2GHz', nchannels=8, stokes='XX YY');


###define the feed
sm.setfeed(mode='perfect X Y',pol=['']);
###
### Limits
sm.setlimits(shadowlimit=0.01, elevationlimit='10deg');
###
### Set autocorrelation weight
sm.setauto(0.0)

###define sources
sm.setfield('mysource', dir0)
sm.setfield('mycal', dircal, 'C')

###now lets observe the cal
time0=me.epoch('TAI','2012/05/21/05:05:13')
sm.settimes(integrationtime='10s', usehourangle=False, 
             referencetime=time0);
##observe cal for 100s ~1 hour before ref time
sm.observe(sourcename='mycal', spwname='band3',
                starttime='-3700s',
                stoptime='-3600s');

### observe source for 2 hours
sm.observe(sourcename='mysource', spwname='band3',
                starttime='-3600s',
                stoptime='3600s');

### observe cal again at the end
sm.observe(sourcename='mycal', spwname='band3',
                starttime='3600s',
                stoptime='3700s');


###now lets put some data in the observation...model of source
sm.setdata(fieldid=0)
sm.predict(imagename=['glmtest.modim'])
###predict cal visibilities
sm.setdata(fieldid=1)
sm.predict(complist='mycal.cl')
sm.done()