u

src/yapsut/__init__.py

@@ -11,6 +11,7 @@ from .profiling import timer,timerDict
 from .struct import struct
 from .DummyCLI import DummyCLI
 from .nearly_even_split import nearly_even_split
+from .geometry import ElipsoidTriangulation 
 
 #from import .grep
 #import .intervalls

src/yapsut/geometry.py

+import numpy as np
+class ElipsoidTriangulation :
+   """ creates an elipsoid triangulated """
+   @property
+   def axis(self) :
+        """ elipsoid axis"""
+        return self._axis
+   @property
+   def shape(self) :
+        """ grid shape"""
+        return self._shape
+   @property
+   def nphi(self) :
+        """ nphi shape"""
+        return self._nphi
+   @property
+   def ntheta(self) :
+        """ ntheta shape"""
+        return self._ntheta
+   @property
+   def phi(self) :
+        """ phi list"""
+        return self._phi
+   @property
+   def theta(self) :
+        """ theta list"""
+        return self._theta
+   @property
+   def simplices(self) :
+        """ list of triangles (simplices)"""
+        return self._Triangles
+   @property
+   def simplices(self) :
+        """ list of triangles (simplices)"""
+        return self._Triangles
+   @property
+   def polarS(self) :
+        """ list of polar coordinates (for each point in each simplice)"""
+        return self._TrAng
+   @property
+   def polarS(self) :
+        """ list of polar coordinates (simplices)"""
+        return self._TrAng
+   @property
+   def vecS(self) :
+        """ list of cartesian coordinates for each triangle"""
+        return self._TrVec
+   @property
+   def verS(self) :
+        """ list of cartesian coordinates for each triangle as versors"""
+        return self._TrVec
+   @property
+   def normalS(self) :
+        """ list of oriented normals for each triangle"""
+        return self._TrNorm
+   @property
+   def polar_normalS(self) :
+        """ list of intercepts for each triangle, polar coordinates"""
+        return self._TrNormAng
+   @property
+   def interceptS(self) :
+        """ list of intercepts for each triangle"""
+        return self._TrItcp
+   @property
+   def areaS(self) :
+        """ list of areas for each triangle"""
+        return self._TrArea
+   @property
+   def areaS(self) :
+        """ list of areas for each triangle"""
+        return self._TrArea
+   @property
+   def vecCpsS(self) :
+        """ list of coplanar coordinates for each triangle"""
+        return self._TrVecCpl
+   #   
+   def __len__(self) :
+        return self._N
+   #   
+   def __init__(self,a,b,c,nphi,ntheta) :
+        self._axis=np.array([a,b,c])
+        self._nphi=nphi
+        self._ntheta=ntheta
+        self._shape=(ntheta,nphi)
+        
+        self._phi=np.linspace(0,1,self._nphi)*2*np.pi
+        self._theta=np.linspace(0,1,self._ntheta)*np.pi
+        
+        phi=self._phi
+        theta=self._theta
+
+        #
+        # map of triangles
+        Triang=[]
+
+        #    north pole triangles
+        itheta=0
+        iP=[0,0]
+        for iiph in np.arange(nphi-1) :
+            iph0=np.mod(iiph,nphi)
+            iph1=np.mod(iiph+1,nphi)
+            #
+            T=[iP,[iph0,itheta+1],[iph1,itheta+1],'0,1']
+            Triang.append(T)
+
+        #    bands triangles
+        Theta=[]
+        for itheta in np.arange(1,ntheta-2) :
+            itheta1=itheta+1
+            ss=str(itheta)+'-'+str(itheta+1)
+            for iiph in np.arange(nphi-1) :
+                iph0=np.mod(iiph,nphi)
+                iph1=np.mod(iiph+1,nphi)
+                #
+                T=[[iph0,itheta],[iph0,itheta1],[iph1,itheta1],ss]
+                Triang.append(T)
+                #
+                T=[[iph1,itheta1],[iph1,itheta],[iph0,itheta],ss]
+                Triang.append(T)
+
+        #    south pole triangles
+        itheta=ntheta-2
+        itheta1=ntheta-1
+        iP=[0,itheta1]
+        ss='0-'+str(itheta1)
+        for iiph in np.arange(nphi-1) :
+            iph0=np.mod(iiph,nphi)
+            iph1=np.mod(iiph+1,nphi)
+            #
+            T=[[iph0,itheta],iP,[iph1,itheta],ss]
+            Triang.append(T)
+        self._Triangles=Triang
+        self._N=len(Triang)
+
+        # convert map of triangles in angles
+        TrAng=[]
+        for BL in Triang :
+            T=[]
+            for iT in BL[:3] : 
+                T.append([phi[iT[0]],theta[iT[1]]])
+            TrAng.append(T)
+        TrAng=np.array(TrAng)
+        self._Angles=TrAng
+
+        # convert map of triangles in vectors
+        TrVec=[]
+        for BL in Triang :
+            T=[]
+            for iT in BL[:3] : 
+                ph=phi[iT[0]]
+                cph=np.cos(ph) ; sph=np.sin(ph) 
+                th=theta[iT[1]]
+                cth=np.cos(th) ; sth=np.sin(th) 
+                T.append([cph*sth,sph*sth,c*cth])
+            TrVec.append(T)
+        TrVec=np.array(TrVec)
+        self._TrVec=TrVec
+
+        # convert map of triangles in versors
+        TrVers=[]
+        for BL in TrVec :
+            T=[]
+            for iT in BL[:3] : 
+                T.append(iT/(iT.dot(iT))**0.5)
+            TrVers.append(T)
+        TrVers=np.array(TrVers)
+        self._TrVers=TrVers
+
+        # normals of triangles as vectors, angles and intercept of planes
+        # TrDprod is for each triangle the n.dot(r_i) with r_i are the triangles vertex
+        TrDprod=[]
+        TrNorm=[]
+        TrNormAng=[]
+        TrItcp=[]
+        for BL in TrVec :
+            r1,r2,r3=BL
+            d21=r2-r1
+            d31=r3-r1
+            cdprod=np.cross(d21,d31)
+            n=cdprod/(cdprod.dot(cdprod))**0.5
+            #
+            sgn=1
+            a1=n.dot(r1)/r1.dot(r1)**0.5
+            a2=n.dot(r2)/r2.dot(r2)**0.5
+            a3=n.dot(r3)/r3.dot(r3)**0.5
+            if a1<0 or a2<0 or a3<0 :
+                n=-1*n
+                sgn=-1
+                a1=n.dot(r1)/r1.dot(r1)**0.5
+                a2=n.dot(r2)/r2.dot(r2)**0.5
+                a3=n.dot(r3)/r3.dot(r3)**0.5
+            TrItcp.append([-n.dot(r1),sgn])
+            TrNorm.append(n)
+            TrDprod.append([a1,a2,a3])
+            #
+            TrNormAng.append([np.arctan2(n[1],n[0]),np.arccos(n[2])])
+        TrItcp=np.array(TrItcp)
+        TrNorm=np.array(TrNorm)
+        TrDprod=np.array(TrDprod)
+        TrNormAng=np.array(TrNormAng)
+        self._TrItcp=TrItcp
+        self._TrNorm=TrNorm
+        self._TrDprod=TrDprod
+        self._TrNormAng=TrNormAng
+
+        # computes the are of triangles and VecCpl, triangles vertex in coplanar representation
+        TrVecCpl=[]
+        TrArea = []
+        for ibl in np.arange(len(TrNorm)) :
+            n=TrNorm[ibl]
+            trv=TrVec[ibl]
+            #
+            nlo,nla=TrNormAng[ibl]
+            cnlo=np.cos(nlo)
+            cnla=np.cos(nla)
+            snlo=np.sin(nlo)
+            snla=np.sin(nla)
+            #
+            Mlo=np.array([[cnlo,-snlo,0],[snlo,cnlo,0],[0,0,1]])
+            Mla=np.array([[cnla,0,snla],[0,1,0],[-snla,0,cnla]])
+            #Mlo.dot(Mla.dot(np.array([0,0,1])))
+            Mlo.dot(Mla.dot(np.array([0,0,1])))
+            nrot=Mla.T.dot(Mlo.T.dot(n))
+            #
+            # method of determinant and closed form
+            # rotate triangle in coplanar coordinates (z = constant)
+            ctr=np.array([Mla.T.dot(Mlo.T.dot(k)) for k in trv])
+            cr1,cr2,cr3=ctr
+            TrVecCpl.append(ctr)
+            # computes area closed form
+            area1=0.5*abs(cr1[0]*(cr2[1]-cr3[1])+cr2[0]*(cr3[1]-cr1[1])+cr3[0]*(cr1[1]-cr2[1]))
+            # computes area method of determinant    
+            ## z forced to be 1
+            #ctr[:,2]=1
+            #area2=0.5*scipy.linalg.det(ctr)
+            #TrArea.append([area2,area1])
+            TrArea.append(area1)
+        TrArea=np.array(TrArea).T
+        TrVecCpl=np.array(TrVecCpl)
+        self._TrArea=TrArea
+        self._TrVecCpl=TrVecCpl
+   #
+   def dotNormals(self,v) :
+      """ computes dot product between normals of simplices and a vector v """
+      return np.array([v.dot(k) for k in self.normalS])
+   #
+   # da verificare
+   #def rotateSimplices(self,lo,la) :
+      #""" computes rotations of simplices  
+          #lo = longitude [rad] 
+          #la = latitude [rad]
+      #"""
+      #cnlo=np.cos(lo)
+      #cnla=np.cos(la)
+      #snlo=np.sin(lo)
+      #snla=np.sin(la)
+      ##
+      #Mlo=np.array([[cnlo,-snlo,0],[snlo,cnlo,0],[0,0,1]])
+      #Mla=np.array([[cnla,0,snla],[0,1,0],[-snla,0,cnla]])
+      #out=[]
+      #for T in self.vecS :
+         #out.append(Mla.T.dot(Mlo.T.dot(k)) for k in T)
+      #return np.array(out)
+   #
+   # da verificare
+   #def rotateNormals(self,lo,la) :
+      #""" computes rotations of normals  
+          #lo = longitude [rad]
+          #la = latitude [rad]
+      #"""
+      #cnlo=np.cos(lo)
+      #cnla=np.cos(la)
+      #snlo=np.sin(lo)
+      #snla=np.sin(la)
+      ##
+      #Mlo=np.array([[cnlo,-snlo,0],[snlo,cnlo,0],[0,0,1]])
+      #Mla=np.array([[cnla,0,snla],[0,1,0],[-snla,0,cnla]])
+      ##
+      #return np.array([Mla.T.dot(Mlo.T.dot(k)) for k in self.normalS])
+