eff_osmose.py
3.91 KB
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
title = "Osmose"
tip = "applies parallactic angle rotation to U and Q"
onein = False
import numpy as np
import nodastro
RAD = np.pi/180
from guidata.qt.QtGui import QMessageBox
from guidata.dataset.datatypes import DataSet
from guidata.dataset.dataitems import (IntItem, StringItem, ChoiceItem, FloatItem, BoolItem)
from guiqwt.config import _
class NOD3_App:
def __init__(self, parent):
self.parent = parent
self.parent.activateWindow()
def compute_app(self, **args):
class FuncParam(DataSet):
FacU = IntItem('FacU', default=1, max=1, min=-1)
FacQ = IntItem('FacQ', default=1, max=1, min=-1)
IAngle = FloatItem('IAngle', default=0.0)
Aver = BoolItem("Average I1 and I2", default=False)
name = title.replace(" ", "")
if args == {}:
param = FuncParam(_("Osmose"), "Reads in U and Q map and corrects for parallactic angle rotation.\nFacU: -1 or 1\nFacQ: -1 or 1\nIAngle: Instrumental Angle")
else:
param = self.parent.ScriptParameter(name, args)
# if no parameter needed set param to None. activate next line
self.parent.compute_11(name, lambda m, p: self.function(m, p), param, onein)
def Error(self, msg):
QMessageBox.critical(self.parent.parent(), title,
_(u"Error:")+"\n%s" % str(msg))
def function(self, ms, p):
if len(ms) < 4:
self.Error("I1, I2, U* and Q* maps are needed")
return [], None
elif len(ms) > 4:
self.Error("Exactly one set of I1, I2, U* and Q* maps are needed")
return [], None
I1 = None
I2 = None
U = None
Q = None
for m in ms:
if not "MAPTYPE" in m.header:
#if not "MAPTYPE" in m.header:
self.Error("MAPTYPE (I,U,Q) not specified in header")
return [], None
else:
if not m.header["MAPTYPE"][0] in ('I', 'U', 'Q'):
self.Error("MAPTYPE (I,U,Q) not specified in header")
return [], None
if m.header["MAPTYPE"] == "I1": I1 = m
if m.header["MAPTYPE"] == "I2": I2 = m
if m.header["MAPTYPE"] == "U": U = m
if m.header["MAPTYPE"] == "Q": Q = m
if type(I1) == type(None) or type(I2) == type(None) or type(U) == type(None) or type(Q) == type(None) :
self.Error("one ore more of MAPTYPE (I,U,Q) are not specified in header")
return [], None
def parmat(p):
return np.array([[np.cos(p), -np.sin(p)], [np.sin(p), np.cos(p)]])
ly, lx = U.parmap.shape
na = nodastro.nodtrafo(m.header['DATE_OBS'])
if 'CROTA2' in m.header: rot = m.header['CROTA2']
#if 'CROTA2' in m.header: rot = m.header['CROTA2']
else: rot = 0.0
if 'EPOCH' in m.header: equinox = m.header['EPOCH']
#if 'EPOCH' in m.header: equinox = m.header['EPOCH']
else: equinox = None
coord = m.header['CTYPE1']
for j in range(ly):
b = m.header['CRVAL2'] + (j-m.header['CRPIX2'])*m.header['CDELT2']
for i in range(lx):
l = m.header['CRVAL1'] + (i-m.header['CRPIX1'])*m.header['CDELT1']
omega = na.getPosAng(l, b, rot, coord, equinox)
eta = (p.IAngle -U.parmap[j][i] - omega)*RAD*2
U.data[j][i], Q.data[j][i] = np.dot(parmat(eta),
[p.FacU*U.data[j][i], p.FacQ*Q.data[j][i]])
if p.Aver:
I1.data = (I1.data + I2.data)/2
I1.header['MAPTYPE'] = 'I'
if "EXTNAME" in I1.header:
#if "EXTNAME" in I1.header:
extname = I1.header["EXTNAME"]
extname = extname.replace("1", "")
I1.header["EXTNAME"] = (extname, "Name of image extension")
return [I1, U, Q], p
else:
return [I1, I2, U, Q], p