KAI doesn't support NIRC2 Vertical Angle Mode in combined images

[nsabrams]

KAI does not currently support vertical angle mode since the position angle is found differently (i.e. from the parallactic angle). So it does not pick up that the images are rotated with respect to one another. For future reference, via Pranav Nagarajan and pyKLIP this is how we potentially could try an incorporate it in the future. (Note pyKLIP apparently has different sign conventions that should be paid attention to).

def get_pa(hdulist, obsdate=None, rotmode=None, mean_PA=True, write_hdr=True):
    """
    Given a FITS data-header unit list (HDUList), returns the NIRC2 PA in [radians].
    PA is angle of detector relative to sky; ROTMODE is rotator tracking mode;
    PARANG is parallactic angle astrometric; INSTANGL is instrument angle;
    ROTPOSN is rotator physical position.
    Additional PA offset of -0.252 or -0.262 deg is applied for NIRC2 narrow cam
    depending on observation date.
    NOTE that the PA sign is flipped at the very end before output to conform to
    pyKLIP's rotation convention.
    
    Inputs:
        hdulist: a FITS HDUList (NOT a single HDU).
        obsdate: date of observation; will try to get from prihdr if not provided.
        rotmode: 'vertical angle' for ADI mode with PA rotating on detector, or
                 'position angle' for mode with PA orientation fixed on detector.
        mean_pa: if True (default), return the mean PA during the exposure.
                 If False, return the PA at the start of the exposure only.
                 Only applies to vertical angle mode.
        write_hdr: if True (default), writes keys to file header and saves them.
    """
    
    prihdr = hdulist[0].header
    
    # Date of NIRC2 servicing that changed PA offset.
    date_2015_4_13 = time.strptime("2015-4-13", "%Y-%m-%d")
    
    # If don't have it, get observation date (UT) from header and make a time object.
    if obsdate is None:
        obsdate = time.strptime(prihdr['DATE-OBS'], "%Y-%m-%d")
    
    # Additional offset to narrow cam PA not included in INSTANGL keyword.
    # This offset changed after instrument servicing on April 13, 2015.
    if prihdr['CAMNAME'].lower() == 'narrow':
        if obsdate < date_2015_4_13:
            zp_offset = -0.252 # [deg]; from Yelda et al. 2010
        elif obsdate >= date_2015_4_13:
            zp_offset = -0.262 # [deg]; from Service et al. 2016
    else:
        zp_offset = 0.
        print("WARNING: No PA offset applied.")
    
    if rotmode is None:
        global _last_rotmode
        try:
            rotmode = prihdr['ROTMODE'].strip().lower()
            _last_rotmode = rotmode
        except:
            rotmode = _last_rotmode.lower()
    rotposn = prihdr['ROTPOSN'] # [deg]
    instangl = prihdr['INSTANGL'] # [deg]

    if rotmode.lower() == 'vertical angle':
        parang = prihdr['PARANG'] # [deg]
        pa_deg = parang + rotposn - instangl + zp_offset # [deg]
    elif rotmode.lower() == 'position angle':
        pa_deg = rotposn - instangl + zp_offset # [deg]
    else:
        raise NotImplementedError
    
    if mean_PA and (rotmode.lower() == 'vertical angle'):
        # Get info for PA smearing calculation.
        epochobj = prihdr['DATE-OBS']
        name = prihdr['targname']
        expref = prihdr['itime']
        coaddref = prihdr['coadds']
        sampref =  prihdr['sampmode']
        msrref =   prihdr['MULTISAM']
        xdimref =  prihdr['naxis1']
        ydimref =  prihdr['naxis2']
        tel = prihdr['TELESCOP']
        dec = Angle(prihdr['DEC'], unit = u.deg).deg + prihdr['DECOFF']
        if tel.lower() == 'keck ii':
            tel = 'keck2' # just cleaning up str
        
        # Calculate total time of exposure (integration + readout).
        if sampref == 2: totexp = ( expref + 0.18*(xdimref/1024.)**2) * coaddref
        if sampref == 3: totexp = ( expref + (msrref-1)*0.18*(xdimref/1024.)**2)*coaddref
        
        tinteg = totexp # [seconds]
        totexp = totexp/3600. # [hours]
        
        # Get hour angle at start of exposure.
        tmpahinit = Angle(prihdr['HA'], unit = u.deg).deg # [deg]
        ahobs = 24.*tmpahinit/360. # [hours]
        
        # Estimate vertical position angle at each second of the exposure.
        vp = [] #fltarr(round(3600.*totexp))
        for j in range(0, int(round(3600.*totexp))-1):
            ahtmp = ahobs + (j*1.+0.001)/3600. # [hours]
            vp.append(par_angle(ahtmp, dec, NIRC2Data.observatory_latitude))
            if j == 0: vpref = vp[0]
        vp = np.array(vp)
        
        # Handle case where PA crosses 0 <--> 360.
        vp[vp < 0.] += 360.
        vp[vp > 360.] -= 360.
        if vpref < 0.: vpref += 360.
        if vpref > 360.: vpref -= 360.
        
        # Check that images near PA=0 are handled correctly.
        if any(vp > 350) & any(vp < 10):
            vp[vp > 350] -= 360
        
        vpmean = np.nanmean(vp)
        
        if (vpmean < 0) & (vpref > 350):
            vpmean += 360.
        
        pa_deg_mean = pa_deg + (vpmean - vpref)
    else:
        pa_deg_mean = pa_deg
        vpmean = np.nan
        vpref = np.nan
        totexp = np.nan
    
    if write_hdr:
        prihdr['TOTEXP'] = (totexp, 'Total exposure time [hours]')
        prihdr['PASTART'] = (pa_deg, "Position angle at exposure start [deg]")
        prihdr['PASMEAR'] = ((vpmean - vpref), "Exposure's weighted-mean PA minus PASTART [deg]")
        prihdr['ROTNORTH'] = (pa_deg_mean, "Mean PA of North during exposure [deg]")
        
        hdulist.flush(output_verify='ignore')
    
    # Flip sign to conform to pyKLIP rotation convention.
    return pa_deg_mean