Commit 98de77cd authored by Michele Maris's avatar Michele Maris
Browse files

u

parent 5e019e20
Loading
Loading
Loading
Loading
+1 −1
Original line number Diff line number Diff line
*.pyc
*__python__*
*__pycache__*
+1 −0
Original line number Diff line number Diff line
2021 dec 01 - Repo created
2021 dec 02 - Uploaded PyPRM 0.0 
2021 dec 03 - Added example
+1 −1
Original line number Diff line number Diff line
@@ -3,4 +3,4 @@ PRM python version
M.Maris
2021 Oct 23
2021 Nov 20
2021 Dec 2
+111 −268
Original line number Diff line number Diff line
"""

__VERSION__="0.0.3 - 2021 Dec 02"
__DESCRIPTION__="""
M.Maris

first version of order 0 Pointing Reconstruction Model (PRM)
@@ -51,206 +53,12 @@ with c=cos(alpha), s=sin(alpha)
Issue 
0.0.0 - 2021 Oct 23 : First implementation
0.0.1 - 2021 Nov 20 : Definition of matrices in agreement with STRIP PRM Doc
0.0.2 - 2021 Nov 22 : Improved defintion of matrices

0.0.2 - 2021 Nov 22 : Improved definition of matrices
0.0.3 - 2021 Dec  2 : Structure improved, added example
"""

import numpy as np

class listMatrix3x3 :
   @property
   def M(self) :
      return self._m
   @M.setter
   def M(self,this) :
      self._m=this
      self._dtype=this.dtype
   @property
   def shape(self) :
      return self._m.shape
   @shape.setter
   def shape(self,this) :
      self._m.shape=this
   @property
   def size(self) :
      return self._m.size
   @property
   def isSingle(self) :
      return len(self.size)==3
   #
   def __init__(self,*arg,dtype=float,value=0) :
      self._dtype=dtype
      if len(arg)==0 :
         self._m=np.array([])
      else :
         if np.isscalar(arg[0]) :
            self._m=self.zeros(arg[0],dtype)+value
         else :
            self._m=np.array(arg[0])
   #
   def __len__(self) :
      return self._m.shape[0]
   #
   def __getitem__(self,this) :
      return self._m[this]
   #
   def __setitem__(self,this,that) :
      self._m[this]=that
   #
   def zeros(self,size,dtype) :
      return np.zeros([size,3,3],dtype=dtype)
   #
   def __clear__(self) :
      self._m=None
   #
   def __str__(self) :
      return self._m.__str__()
   #
   def dot(self,this) :
      if np.isscalar(this) :
         out=self.__class__(self._m.dot(this))
         return out
      out=[]
      for k2 in self._m :
         for k1 in this :
            out.append(k2.dot(k1))
      return self.__class__(out)
   #
   def apply_to(self,this_lm) :
      if len(this_lm)==len(self) :
         out=listMatrix3x3(len(self))
         for ik in range(len(self)) :
            out[ik]=self[ik].dot(this_lm[ik])
         return out
      elif len(this_lm)==1 and len(self)>1 :
         out=listMatrix3x3(len(this_lm))
         for ik in range(len(this_lm)) :
            out[ik]=self[0].dot(this_lm[ik])
         return out
      elif len(this_lm)>1 and len(self)==1 :
         out=listMatrix3x3(len(self))
         for ik in range(len(self)) :
            out[ik]=self[ik].dot(this_lm[0])
         return out
      else :
         raise Exception("error: uncompatible dimesions between arrays","")
   ##
   #def __mul__(self,this) :
      #if np.isscalar(this) :
         #out=self.__class__(this*self._m)
         #return out
      ##
      ##if self.__class__==this.__class__:
         ##out=[]
         ##for k2 in self._m :
            ##for k1 in this.M :
               ##out.append(k2.dot(k1))
         ##return self.__class__(out)
      ###
      ##ooo=str(self.__class__).split('.')[-1].split("'")[0]
      ##raise Exception("error: a "+ooo+" object can be multiplied by a scalar or another "+ooo+" object ","")
      ##
      ##if this.__class__==self._m.__class__:
         ##print("__mul__ 1")
         ##out=self.__class__([k.dot(this) for k in self._m])
         ##return out
   ##
   #def __rmul__(self,this) :
      #if np.isscalar(this) :
         #return self.__class__(this*self._m)
      ###
      ##if self.__class__==this.__class__:
         ##out=[]
         ##for k1 in this.M :
            ##for k2 in self._m :
               ##out.append(k2.dot(k1))
         ##return self.__class__(out)
      ###
      ##ooo=str(self.__class__).split('.')[-1].split("'")[0]
      ##raise Exception("error: a "+ooo+" object can be multiplied by a scalar or another "+ooo+" object ","")
      ##
      ##if this.__class__==self._m.__class__:
         ##out=self.__class__([this.dot(k) for k in self._m])
         ##return out
      #
   def copy(self) :
      import copy
      return copy.deepcopy(self)
   def slice(self,*arg) :
      i1=arg[0]
      if len(arg)==2 :
         i2=arg[1]
      else :
         i2=i1+1
      out=listMatrix3x3()
      out._m=self._m[i1:i2]
      return out

#m=np.array([[1,0,0],[1,1,0],[1,1,1]])
#lm1=listMatrix3x3([m])
#lm=listMatrix3x3([k*m for k in range(1,100)])

class ElemRot(listMatrix3x3) :
   @property
   def axis(self) :
      return self._axis
   @property
   def angle(self) :
      return self._angle
   @angle.setter
   def angle(self,this) :
      self._angle=np.array([this]) if np.isscalar(this) else this
      self._c=np.cos(self._angle)
      self._s=np.sin(self._angle)
      self._m=self.zeros(len(self._angle),float)
      self._m[:,0,0]=1
      self._m[:,1,1]=1
      self._m[:,2,2]=1
      #
      if self.axis=='x' :
         self._m[:,1,1]=self._c
         self._m[:,1,2]=-self._s
         self._m[:,2,1]=self._s
         self._m[:,2,2]=self._c
      elif self.axis=='y' :
         self._m[:,0,0]=self._c
         self._m[:,0,2]=self._s
         self._m[:,2,0]=-self._s
         self._m[:,2,2]=self._c
      elif self.axis=='z' :
         self._m[:,0,0]=self._c
         self._m[:,0,1]=-self._s
         self._m[:,1,0]=self._s
         self._m[:,1,1]=self._c
      else :
         raise Exception("error: undefined rotation axis","")
   def __init__(self,axis,Angle,deg=True) :
      listMatrix3x3.__init__(self)
      self._axis=None
      self._angle=None
      if axis==0 or str(axis).lower()=='x' :
         self._axis='x'
      elif axis==1 or str(axis).lower()=='y' :
         self._axis='y'
      elif axis==2 or str(axis).lower()=='z' :
         self._axis='z'
      elif type(axis)==type("") :
         self._axis=str(axis)
      else :
         raise Exception("error: undefined rotation axis","")
      if not Angle is None :
         if deg :
            self.angle=np.deg2rad(Angle)
         else :
            self.angle=Angle

#erot=ElemRot('x',None)
#erot1=ElemRot('x',np.arange(0,10))
#m=np.array([[1,0,0],[1,1,0],[1,1,1]])
#lm1=listMatrix3x3([m])
#lm=listMatrix3x3([k*m for k in range(1,100)])
#stop

class PRM0 :
   @property
   def roll(self) :
@@ -259,7 +67,7 @@ class PRM0 :
   def roll(self,this) :
      self._d['roll']=this
      self._mark['roll']=1
      self._matr['roll']=self.MatrZ(np.deg2rad(np.array([this])))
      self._matr['roll']=self.ERotZ(np.deg2rad(np.array([this])))
      self._isComplete=False
   #
   @property
@@ -269,7 +77,7 @@ class PRM0 :
   def pan(self,this) :
      self._d['pan']=this
      self._mark['pan']=1
      self._matr['pan']=self.MatrY(np.deg2rad(np.array([this])))
      self._matr['pan']=self.ERotY(np.deg2rad(np.array([this])))
      self._isComplete=False
   #
   @property
@@ -279,7 +87,7 @@ class PRM0 :
   def tilt(self,this) :
      self._d['tilt']=this
      self._mark['tilt']=1
      self._matr['tilt']=self.MatrX(np.deg2rad(np.array([this])))
      self._matr['tilt']=self.ERotX(np.deg2rad(np.array([this])))
      self._isComplete=False
   #
   @property
@@ -289,7 +97,7 @@ class PRM0 :
   def ctheta0(self,this) :
      self._d['ctheta0']=this
      self._mark['ctheta0']=1
      self._matr['ctheta0']=self.MatrY(-np.deg2rad(np.array([this])))
      self._matr['ctheta0']=self.ERotY(-np.deg2rad(np.array([this])))
      self._isComplete=False
   #
   @property
@@ -299,7 +107,7 @@ class PRM0 :
   def cphi0(self,this) :
      self._d['cphi0']=this
      self._mark['cphi0']=1
      self._matr['cphi0']=self.MatrZ(-np.deg2rad(np.array([this])))
      self._matr['cphi0']=self.ERotZ(-np.deg2rad(np.array([this])))
      self._isComplete=False
   #
   @property
@@ -309,27 +117,27 @@ class PRM0 :
   def tiltFork(self,this) :
      self._d['tiltFork']=this
      self._mark['tiltFork']=1
      self._matr['tiltFork']=self.MatrX(np.deg2rad(np.array([this])))
      self._matr['tiltFork']=self.ERotX(np.deg2rad(np.array([this])))
      self._isComplete=False
   #
   @property
   def zVax(self) :
      return self._d['zVax']
   @zVax.setter
   def zVax(self,this) :
      self._d['zVax']=this
      self._mark['zVax']=1
      self._matr['zVax']=self.MatrX(np.deg2rad(np.array([this])))
   def zVAX(self) :
      return self._d['zVAX']
   @zVAX.setter
   def zVAX(self,this) :
      self._d['zVAX']=this
      self._mark['zVAX']=1
      self._matr['zVAX']=self.ERotX(np.deg2rad(np.array([this])))
      self._isComplete=False
   #
   @property
   def omegaVax(self) :
      return self._d['omegaVax']
   @omegaVax.setter
   def omegaVax(self,this) :
      self._d['omegaVax']=this
      self._mark['omegaVax']=1
      self._matr['omegaVax']=self.MatrZ(np.deg2rad(np.array([this])))
   def omegaVAX(self) :
      return self._d['omegaVAX']
   @omegaVAX.setter
   def omegaVAX(self,this) :
      self._d['omegaVAX']=this
      self._mark['omegaVAX']=1
      self._matr['omegaVAX']=self.ERotZ(np.deg2rad(np.array([this])))
      self._isComplete=False
   #
   @property
@@ -343,19 +151,6 @@ class PRM0 :
            return False
      return True
   #
   def keys(self) :
      return ['roll','pan','tilt','ctheta0','cphi0','tiltFork','zVax','omegaVax']
   #
   def __getitem__(self,this) :
      return self._d[this]
   #
   def __setitem__(self,this,that) :
      self._d[this]=that
   #
   def copy(self) :
      import copy
      return copy.deepcopy(self)
   #
   def __init__(self,ideal=False) :
      self._isComplete=False
      self._d={}
@@ -370,30 +165,88 @@ class PRM0 :
      if ideal :
         self.setIdeal()
   #
   def keys(self) :
      return ['roll','pan','tilt','ctheta0','cphi0','tiltFork','zVAX','omegaVAX']
   #
   def __getitem__(self,this) :
      return self._d[this]
   #
   def __setitem__(self,this,that) :
      self._d[this]=that
   #
   def M(self,k) :
      return self._matr[k]
   #
   def copy(self) :
      import copy
      return copy.deepcopy(self)
   #
   def _base_matrix(self,Angle) :
      _angle=np.array([Angle]) if np.isscalar(Angle) else Angle
      c=np.cos(_angle)
      s=np.sin(_angle)
      I=np.zeros([len(_angle),3,3])
      I[:,0,0]=1
      I[:,1,1]=1
      I[:,2,2]=1
      return I,c,s
   #
   def ERotX(self,Angle) :
      R,c,s=self._base_matrix(Angle)
      R[:,1,1]= c
      R[:,1,2]=-s
      R[:,2,1]= s
      R[:,2,2]= c
      return R
   #
   def ERotY(self,Angle) :
      R,c,s=self._base_matrix(Angle)
      R[:,0,0]= c
      R[:,0,2]= s
      R[:,2,0]=-s
      R[:,2,2]= c
      return R
   #
   def ERotZ(self,Angle) :
      R,c,s=self._base_matrix(Angle)
      R[:,0,0]= c
      R[:,0,1]=-s
      R[:,1,0]= s
      R[:,1,1]= c
      return R
   #
   def setIdeal(self) :
      self.roll=0.
      self.pan=0.
      self.tilt=0.
      self.tiltFork=0.
      self.zVax=0.
      self.omegaVax=0.
      self.zVAX=0.
      self.omegaVAX=0.
      self.cphi0=0.
      self.ctheta0=0.
      self.commit()
   #
   def commit(self) :
      if self.isAllSet :
         out=self._matr['pan'].apply_to(self._matr['roll'])
         out=self._matr['tilt'].apply_to(out)
         self._matr['M1']=self._matr['ctheta0'].apply_to(out)
         self._matr['M2']=self._matr['cphi0'].apply_to(self._matr['tiltFork'])
         self._matr['M3']=self._matr['omegaVax'].apply_to(self._matr['zVax'])
         out=self._matr['tilt'][0].dot(self._matr['pan'][0].dot(self._matr['roll'][0]))
         out=self._matr['ctheta0'][0].dot(out)
         out.shape=(1,3,3)
         self._matr['M1']=out*1
         #
         out=self._matr['cphi0'][0].dot(self._matr['tiltFork'][0])
         out.shape=(1,3,3)
         self._matr['M2']=out*1
         #
         out=self._matr['omegaVAX'][0].dot(self._matr['zVAX'][0])
         out.shape=(1,3,3)
         self._matr['M3']=out*1
         #
         self._isComplete=True
      else :
         self._isComplete=False
         raise Error('Not configuration not completed','')
         raise Error('Configuration not completed','')
   #
   def __call__(self,CTHETA_PHI,axis0isNPOS=False) :
   def __call__(self,CTHETA_PHI) :
      """ 
Apply rotations for a list of NPOS control angles: (CTR_THETA, CTR_PHI)

@@ -441,39 +294,29 @@ while PMatrB[0] is the rotation matrix for the case (45.,0.) and PMatrB[:,2] is
      """ 
      if not self.isComplete :
         raise Error('Configuration not completed or not committed','')
      self._d['CTHETA_PHI']=CTHETA_PHI
      #
      self._matr['CTRL_THETA']=self.MatrY(np.deg2rad(CTHETA_PHI[:,0]))
      self._matr['CTRL_PHI']=self.MatrZ(np.deg2rad(CTHETA_PHI[:,1]))
      self._d['CTHETA_PHI']=CTHETA_PHI
      #
      out=self._matr['CTRL_THETA'].apply_to(self._matr['M1'])
      out=self._matr['M2'].apply_to(out)
      out=self._matr['CTRL_PHI'].apply_to(out)
      out=self._matr['M3'].apply_to(out)
      # theta rotation
      self._matr['CTRL_THETA']=self.ERotY(np.deg2rad(CTHETA_PHI[:,0]))
      #
      #stores Rout 
      self._matr['Rout']=out
      # phi rotation
      self._matr['CTRL_PHI']=self.ERotZ(np.deg2rad(CTHETA_PHI[:,1]))
      #
      #set output has NPOS,3,3 shape if required
      #if axis0isNPOS==True :
         #out=np.moveaxis(out,0,2)
      # apply each theta rotation to M1
      out=[k.dot(self._matr['M1'][0]) for k in self._matr['CTRL_THETA']]
      #
      return out
      # apply M2 to each out
      out1=[self._matr['M2'][0].dot(k) for k in out]
      #
   def MatrX(self,Angle,deg=True) :
      return ElemRot('x',Angle,deg=deg)
      # apply phi rotation to each out2
      out2=[self._matr['CTRL_PHI'][ik].dot(k) for ik, k in enumerate(out1)]
      #
   def MatrY(self,angle) :
      return ElemRot('y',Angle,deg=deg)
      # apply M3 to each out3
      out3=np.array([self._matr['M3'][0].dot(k) for k in out2])
      #
   def MatrZ(self,angle) :
      return ElemRot('z',Angle,deg=deg)
      #stores Rout 
      self._matr['Rout']=out3
      #
   def rmult(self,this,that) :
      out=np.zeros([3,3,max(this.shape[2],that.shape[2])])
      for i in range(3) :
         for j in range(3) :
            for l in range(3) :
               out[i,j]+=this[i,l]*that[l,j]
      return out
      return out3
      
+128 −0
Original line number Diff line number Diff line
#
# example of use of pyPRM
#
from PyPRM import *

from collections import OrderedDict
from matplotlib import pyplot as plt

# class to create standardized plots
class stereographic :
   def __init__(self,figsize=[12,12]) :
      self._uxyz=None
      self.h=None
   def project(self,xyz) :
      theta=np.arctan2(xyz[:,1],xyz[:,0])
      radius=np.arccos(xyz[:,2])
      return theta,radius
   @property
   def uxyz(self) :
      return self._uxyz
   @uxyz.setter
   def uxyz(self,this) :
      self._uxyz=this
   def uplot(self) :
      x=np.arange(10)/10.
      self.h=plt.polar(x*2*pi/10,x*0+1.2)
      self.h[0].set_visible(False)
      for k in self.uxyz :
         _utheta,_uradius=self.project(k)
         plt.polar(_utheta,_uradius,color='gray')
   def plot(self,xyz,**karg) :
      theta,radius=self.project(xyz)
      plt.polar(theta[0],radius[0],'o',markersize=8,color='gray')
      plt.polar(theta,radius,**karg)
      plt.show()

# 
# list of Theta,Phi positions
CT_P=np.array([
      [0.,0.]
   ,[30.,0.]
   ,[60.,0.]
   ,[30.,90.]
   ,[30.,180.]
   ,[30.,270.]
   ])

# the first case, an ideal telescope
pm0I=PRM0(ideal=True)
R=OrderedDict()
R['Ideal']=pm0I(CT_P)

# camera roll
pm0roll=PRM0(ideal=True)
pm0roll.roll=45
pm0roll.commit()
R['roll']=pm0roll(CT_P)

# camera pan
pm0pan=PRM0(ideal=True)
pm0pan.pan=10
pm0pan.commit()
R['pan']=pm0pan(CT_P)

# camera tilt
pm0tilt=PRM0(ideal=True)
pm0tilt.tilt=10
pm0tilt.commit()
R['tilt']=pm0tilt(CT_P)

# tiltFork
pm0tiltFork=PRM0(ideal=True)
pm0tiltFork.tiltFork=10
pm0tiltFork.commit()
R['tiltFork']=pm0tiltFork(CT_P)

# ctheta0
pm0ctheta0=PRM0(ideal=True)
pm0ctheta0.ctheta0=10
pm0ctheta0.commit()
R['ctheta0']=pm0ctheta0(CT_P)

# cphi0
pm0cphi0=PRM0(ideal=True)
pm0cphi0.cphi0=10
pm0cphi0.commit()
R['cphi0']=pm0cphi0(CT_P)

# zVAX
pm0zVAX=PRM0(ideal=True)
pm0zVAX.zVAX=10
pm0zVAX.commit()
R['zVAX']=pm0zVAX(CT_P)

# omegaVAX
pm0omegaVAX=PRM0(ideal=True)
pm0omegaVAX.omegaVAX=10
pm0omegaVAX.commit()
R['omegaVAX']=pm0omegaVAX(CT_P)

#
# rotates a 4 points fow
# FOW trace
w=6; sw = np.sin(np.deg2rad(w))
h=2; sh = np.sin(np.deg2rad(h))
z=(1-sw**2-sh**2)**0.5
V=np.array([
    [ sw, sh, z]
   ,[ sw,-sh, z]
   ,[-sw,-sh, z]
   ,[-sw, sh, z]
   ,[ sw, sh, z]
   ])

plt.close('all')
#
std=stereographic()
std.uxyz=[np.array([k.dot(k1) for k1 in V]) for k in R['Ideal']]
#
for karg in ['Ideal','roll','pan','tilt','tiltFork','ctheta0','cphi0','zVAX','omegaVAX'] :
   plt.figure(figsize=[8,8])
   std.uplot()
   for k in R[karg] :
      Vr=np.array([k.dot(k1) for k1 in V])
      std.plot(Vr)
   plt.title(karg)
   plt.savefig('/tmp/polar_%s.pdf'%karg,dpi=plt.gcf().get_dpi())
   plt.show()