CSFPA/CSFPA_main.py at master · seamusomurchu/CSFPA · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
import qubic
import glob
import numpy as np
import os
import timeit
from CSFPA_dataIO import getXYcoords, dataIO, dataAnalysis, SaveVars, IntensityCalc, IntensityCalcRAW

def MainProg(filename):
    start = timeit.default_timer()

    rep = '/home/james/files4CSFPA/Fromqbdataio/'
    repfile = rep + filename
    print "file and rep: ", repfile
    #scene = QubicScene(256)
    #inst = instrument.QubicInstrument(filter_nu=150e9)
#    inst = instrument.QubicInstrument()
#    centers = inst.horn.center[:,0:2]
#    detcenters = inst.detector.center[..., :2]
    d = qubic.qubicdict.qubicDict()
    d.read_from_file('/home/james/qubic/qubic/scripts/global_default_FI.dict')

    q = qubic.QubicMultibandInstrument(d)

#    rep = "/home/james/CF-Source-Focal-Place-Analysis/UItest"
#
#    files = glob.glob(rep+filename)
#    files.sort()

    #vtxs = vertexes of TD detectors
#    vtxs = inst.detector.vertex[496:744]
#    #full detector setup
#    vtxsF = inst.detector.vertex

    #SET DETECTOR VERTEXES TO FULL INSTRUMENT
    #q[0].calibration.detarray = '/usr/local/lib/python2.7/dist-packages/qubic/calfiles/CalQubic_DetArray_FI.fits'
    vtxs = q[0].detector.vertex

    #should be 992 for FI
    #vtxcounter = np.zeros(248)
    vtxcounter = np.zeros(992)
    #DEBUG
    #print "vtxs = ", vtxs.shape
    #print "vtxcounter = ", vtxcounter.shape


    MagXarr, PhaXarr, ReXarr, ImXarr, MagYarr, PhaYarr, ReYarr, ImYarr, vtxcntarr, PixCenX, PixCenY = getXYcoords(repfile,vtxs)
    print "getxycoordfunctest", max(MagXarr), MagXarr.shape

    vtxcounter = np.vstack((vtxcounter,vtxcntarr))
    vtxcounter = vtxcounter.T
    vtxcounter = vtxcounter[:,1:3]

    #IntX = (ReXarr*ReXarr) + (ImXarr*ImXarr)
    #IntY = (ReYarr*ReYarr) + (ImYarr*ImYarr)
    #IntT = IntX[:] + IntY[:]

	#outsourcing intensity calc to dataIO
    IntX, IntY, IntT = IntensityCalc(MagXarr, PhaXarr, MagYarr, PhaYarr)
    print "intensity tests shape max", IntX.shape, max(IntX)
    #use this order for a header
    #dat = np.hstack((MagXmat,vtxcounter,ReXmat,ImXmat,ReYmat,ImYmat))
    dat = np.vstack((MagXarr, PhaXarr, ReXarr, ImXarr, MagYarr, PhaYarr, ReYarr, ImYarr, vtxcntarr, PixCenX, PixCenY, IntX, IntY, IntT))
    dat = dat.T
    #SAVE DATA TO MAKE LOADING EASY
    #create strings for saving to file IO
    datstring = "dat"
    dataIO(dat,datstring)

    datmodstring = "datmod"
    datmod = dataAnalysis(dat) # THIS DIVIDES BY NUMBER OF POINTS PER PIXEL
    dataIO(datmod,datmodstring)

    #load data for plotting raw modal data
    dataCF1 = np.loadtxt(repfile, skiprows=1)
    xycoords = np.array(dataCF1[:,2:4])

    #plotting RAW INTENSITY
    #Ix = (dataCF1[:,4]*np.cos(dataCF1[:,5]))**2 + (dataCF1[:,4]*np.sin(dataCF1[:,5]))**2
    #Iy = (dataCF1[:,6]*np.cos(dataCF1[:,7]))**2 + (dataCF1[:,6]*np.sin(dataCF1[:,7]))**2
    #IT = Ix + Iy

    Ix, Iy, IT = IntensityCalcRAW(repfile)

    #testing setting zeros to NANs
    ITnans = [np.nan if x == 0 else x for x in IT]
    ITnans = np.asarray(ITnans)

    #save vars as local
    SaveVars(MagXarr, PhaXarr, ReXarr, ImXarr, MagYarr, PhaYarr, ReYarr, ImYarr, vtxcntarr, PixCenX, PixCenY, IntX, IntY, IntT, Ix, Iy, IT, xycoords, filename)

    os.system('spd-say "Main program has finished"')
    stop = timeit.default_timer()
    time = stop - start
    seconds = (time - int(time)) * 60
    print time/60, "m", seconds, "s"

    return