Phase binning code for producing pseudomaps

phasebin.py

import numpy as np
import cPickle
import pycfitsio
from collections import defaultdict
import healpy

from planck import Planck

from testenv import remix, todtools
from planck.pointingtools import AHF_btw_OBT

from planck import Planck, pointing, metadata


def unwrap(phase):
    i,=np.where(np.diff(phase)<0)
    i+=1
    offsets = 360. * np.repeat(np.arange(len(i)+1), np.concatenate([i, [len(phase)]]) -np.concatenate([[0],i]))
    return phase + offsets

def do1od(od, freq, cal):
    BINS = {30:1950,
    44:1950,
    70:1950}

    folder = '/global/scratch/sd/planck/user/zonca/data/LFI_%s_hrflag_conv' % cal
    ##### CHANNEL AND OBT
    f = Planck.Planck().f[freq]

    ds = metadata.DataSelector(None)
    PP = ds.get_PP(od)

    phasebinned=defaultdict(dict)
    print(od)
    ahf_files = AHF_btw_OBT([PP[0].start, PP[-1].stop])
    phase = np.array([])
    ahf_obt = np.array([])
    for ahf_file in ahf_files:
        print(ahf_file)
        with pycfitsio.File(ahf_file) as ahf_fitsfile:
            phase = np.concatenate([phase, ahf_fitsfile[0].read_column("SPIN_PHASE")])
            ahf_obt = np.concatenate([ahf_obt, ahf_fitsfile[0].read_column("OBT_SPL")/2.**16])
    phase = unwrap(phase)
    phasebin_edges = np.linspace(0, 360, BINS[freq]+1)
    with pycfitsio.open(folder + '/%04d/L%03d-%04d-R-20110101.fits' % (od, freq, od)) as dipolefile:
        obt = dipolefile['OBT'].read_column('OBT')/2.**16
        baselinefilename = '/project/projectdirs/planck/data/mission/dx7/calibration_run/%s/%d_91-640_1024_60s_baselines.fits' % (cal, freq)
        print('Open %s' % baselinefilename)
        bas = pycfitsio.read(baselinefilename)
        for ch in f.ch:
            print(ch)
            dipole = dipolefile[ch.tag].read_column(ch.tag) - todtools.expand_calibration(bas['START'], bas[ch.tag], obt)
            for PID in PP:
                if PID.stop - PID.start > 60:
                    timerange = np.logical_and(obt>=PID.start, obt<=PID.stop)
                    obt_pp = obt[timerange]
                    dipole_pp = dipole[timerange]
                    assert np.all(np.diff(ahf_obt) >= 0)
                    phase_pp = np.mod(np.interp(obt_pp, ahf_obt, phase), 360.)

                    phasebinned[ch.tag][PID.number] /= hits
                    #print("Empty bins %d" % (hits==0).sum())
                    #dipole[timerange]

    with open('/project/projectdirs/planck/user/zonca/issues/calrun/pmdiff/pkl%d/phasebin_%s_%03d.pkl' % (freq, cal, od), 'w') as f:
        cPickle.dump(phasebinned, f, -1)

if __name__ == '__main__':
    do1od(100, 30, '369S')