i2curv_rep.F File Reference

#include "implicit_f.inc"
#include "com08_c.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	i2curv_rep (inod, x, v, ls1, ls2, lt1, lt2, e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z)

Function/Subroutine Documentation

i2curv_rep()

subroutine i2curv_rep	(	integer, dimension(4)	inod,
			x,
			v,
			ls1,
			ls2,
			lt1,
			lt2,
			e1x,
			e1y,
			e1z,
			e2x,
			e2y,
			e2z,
			e3x,
			e3y,
			e3z )

Definition at line 30 of file i2curv_rep.F.

C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER INOD(4)
C     REAL
      my_real
     .   x(3,*),v(3,*)
      my_real
     .   e1x,e1y,e1z,e2x,e2y,e2z,e3x,e3y,e3z,
     .   ls1,ls2,lt1,lt2
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com08_c.inc"
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER J, JJ, NIR
C     REAL
      my_real
     .   x21,x32,x34,x41,y21,y32,y34,y41,z21,z32,z34,z41,
     . tp22223333,suma,s1,s2
      my_real
     .   xm(4),ym(4),zm(4)
C=======================================================================
      nir=4                                                           
      IF (inod(3) == inod(4)) THEN                                      
        nir=3                                                         
        xm(4)=zero                                                    
        ym(4)=zero                                                    
        zm(4)=zero                                                    
      ENDIF                                                           
      DO jj=1,nir                                                     
        j=inod(jj)  
        xm(jj)=x(1,j)+ v(1,j)*dt12*half                                                
        ym(jj)=x(2,j)+ v(2,j)*dt12*half                                                
        zm(jj)=x(3,j)+ v(3,j)*dt12*half                                                
      ENDDO                                                           
C-----local frame main                                             
      x21 = xm(2) - xm(1)                                             
      x32 = xm(3) - xm(2)                                             
      x34 = xm(3) - xm(4)                                             
      x41 = xm(4) - xm(1)                                             
      y21 = ym(2) - ym(1)                                             
      y32 = ym(3) - ym(2)                                             
      y34 = ym(3) - ym(4)                                             
      y41 = ym(4) - ym(1)                                             
      z21 = zm(2) - zm(1)                                             
      z32 = zm(3) - zm(2)                                             
      z34 = zm(3) - zm(4)                                             
      z41 = zm(4) - zm(1)                                             
                                                                      
      e1x = (x21+x34 )                                                
      e1y = (y21+y34 )                                                
      e1z = (z21+z34 )                                                
 
      e2x = (x32+x41 )                                                
      e2y = (y32+y41 )                                                
      e2z = (z32+z41 )                                                
 
      e3x = e1y*e2z-e1z*e2y                                           
      e3y = e1z*e2x-e1x*e2z                                           
      e3z = e1x*e2y-e1y*e2x                                           
C
      suma   = e3x*e3x+e3y*e3y+e3z*e3z                                
      suma   = one / max(sqrt(suma),em20)                              
      e3x = e3x * suma                                                
      e3y = e3y * suma                                                
      e3z = e3z * suma                                                
C
      s1     = e1x*e1x+e1y*e1y+e1z*e1z                                
      s2     = e2x*e2x+e2y*e2y+e2z*e2z                                
      suma   = sqrt(s1/s2)                                            
      e1x = e1x + (e2y*e3z-e2z*e3y)*suma                              
      e1y = e1y + (e2z*e3x-e2x*e3z)*suma                              
      e1z = e1z + (e2x*e3y-e2y*e3x)*suma                              
C
      suma   = e1x*e1x+e1y*e1y+e1z*e1z                                
      suma   = one / max(sqrt(suma),em20)                              
      e1x = e1x * suma                                                
      e1y = e1y * suma                                                
      e1z = e1z * suma                                                
C
      e2x = e3y * e1z - e3z * e1y                                     
      e2y = e3z * e1x - e3x * e1z                                     
      e2z = e3x * e1y - e3y * e1x                                     
C----- longueurs                                            
      ls1 = sqrt(x21**2 + y21**2 + z21**2)                  
      ls2 = sqrt(x34**2 + y34**2 + z34**2)                  
      lt1 = sqrt(x41**2 + y41**2 + z41**2)                  
      lt2 = sqrt(x32**2 + y32**2 + z32**2)   
C-----------
      RETURN