ov_comp.f

Go to the documentation of this file.

!                   ******************
                    SUBROUTINE ov_comp
!                   ******************
!
     & ( op , x , y , z , c , npoin,
     &   x_err , y_err , z_err)
!
!***********************************************************************
! BIEF   V6P1                                   21/08/2010
!***********************************************************************
!
!brief    OPERATIONS ON VECTORS.
!code
!+   OP IS A STRING OF 8 CHARACTERS, WHICH INDICATES THE OPERATION TO BE
!+   PERFORMED ON VECTORS X,Y AND Z AND CONSTANT C.
!+
!+   THE RESULT IS VECTOR X.
!+
!+   OP = 'X=C     '     :  SETS X TO C
!+   OP = 'X=Y     '     :  COPIES Y IN X
!+   OP = 'X=+Y    '     :  IDEM
!+   OP = 'X=-Y    '     :  COPIES -Y IN X
!+   OP = 'X=1/Y   '     :  COPIES INVERSE OF Y IN X
!+   OP = 'X=Y+Z   '     :  ADDS Y AND Z
!+   OP = 'X=Y-Z   '     :  SUBTRACTS Z TO Y
!+   OP = 'X=YZ    '     :  Y.Z
!+   OP = 'X=-YZ   '     :  -Y.Z
!+   OP = 'X=XY    '     :  X.Y
!+   OP = 'X=X+YZ  '     :  ADDS Y.Z TO X
!+   OP = 'X=X-YZ  '     :  SUBTRACTS Y.Z FROM X
!+   OP = 'X=CXY   '     :  C.X.Y
!+   OP = 'X=CYZ   '     :  C.Y.Z
!+   OP = 'X=CXYZ  '     :  C.X.Y.Z
!+   OP = 'X=X+CYZ '     :  ADDS C.Y.Z TO X
!+   OP = 'X=Y/Z   '     :  DIVIDES Y BY Z
!+   OP = 'X=CY/Z  '     :  DIVIDES C.Y BY Z
!+   OP = 'X=CXY/Z '     :  DIVIDES C.X.Y BY Z
!+   OP = 'X=X+CY/Z'     :  ADDS C.Y/Z TO X
!+   OP = 'X=X+Y   '     :  ADDS Y TO X
!+   OP = 'X=X-Y   '     :  SUBTRACTS Y FROM X
!+   OP = 'X=CX    '     :  MULTIPLIES X BY C
!+   OP = 'X=CY    '     :  MULTIPLIES Y BY C
!+   OP = 'X=Y+CZ  '     :  ADDS C.Z TO Y
!+   OP = 'X=X+CY  '     :  ADDS C.Y TO X
!+   OP = 'X=SQR(Y)'     :  SQUARE ROOT OF Y
!+   OP = 'X=ABS(Y)'     :  ABSOLUTE VALUE OF Y
!+   OP = 'X=N(Y,Z)'     :  NORM OF THE VECTOR WITH COMPONENTS Y AND Z
!+   OP = 'X=Y+C   '     :  ADDS C TO Y
!+   OP = 'X=X+C   '     :  ADDS C TO X
!+   OP = 'X=Y**C  '     :  Y TO THE POWER C
!+   OP = 'X=COS(Y)'     :  COSINE OF Y
!+   OP = 'X=SIN(Y)'     :  SINE OF Y
!+   OP = 'X=ATN(Y)'     :  ARC TANGENT OF Y
!+   OP = 'X=A(Y,Z)'     :  ATAN2(Y,Z)
!+   OP = 'X=+(Y,C)'     :  MAXIMUM OF Y AND C
!+   OP = 'X=-(Y,C)'     :  MINIMUM OF Y AND C
!+   OP = 'X=+(Y,Z)'     :  MAXIMUM OF Y AND Z
!+   OP = 'X=-(Y,Z)'     :  MINIMUM OF Y AND Z
!+   OP = 'X=YIFZ<C'     :  COPIES Y IN X IF Z < C , FOR EACH POINT
!+   OP = 'X=C(Y-Z)'     :  MULTIPLIES C BY (Y-Z)
!
!warning  DIVIDE OPERATIONS INTERNALLY TAKE CARE OF DIVISIONS BY 0.
!+            SUCCESSFUL EXIT OF OV IS THEREFORE NOT A PROOF THAT Y OR Z
!+            NEVER ARE 0. FOR SUCH OPERATIONS, PREFERABLY USE OVD
!
!history  J-M HERVOUET (LNHE)     ; F  LEPEINTRE (LNH)
!+        17/03/2010
!+        V6P0
!+   ORIGINAL IDEA IN ULYSSE. THANK YOU D. LAURENCE
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| C              |-->| A GIVEN CONSTANT
!| NPOIN          |-->| SIZE OF VECTORS
!| OP             |-->| STRING INDICATING THE OPERATION TO BE DONE
!| X              |<--| RESULTING VECTOR
!| Y              |-->| TO BE USED IN THE OPERATION
!| Z              |-->| TO BE USED IN THE OPERATION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)             :: NPOIN
      DOUBLE PRECISION, INTENT(IN)    :: Y(npoin),Z(*),C
      DOUBLE PRECISION, INTENT(INOUT) :: X(npoin)
      DOUBLE PRECISION,OPTIONAL,INTENT(IN):: Y_ERR(npoin),Z_ERR(npoin)
      DOUBLE PRECISION,OPTIONAL, INTENT(INOUT) :: X_ERR(npoin)
      CHARACTER(LEN=8), INTENT(IN)    :: OP
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER I
      DOUBLE PRECISION ERROR,ERROR1,TEMP(npoin),TEST
      DOUBLE PRECISION TMP
!
      INTRINSIC sqrt,abs,cos,sin,atan,max,min
!
!-----------------------------------------------------------------------
!
      temp=0.d0
      error=0.d0
      error1=0.d0
      SELECT CASE(op(3:8))
!
!-----------------------------------------------------------------------
!
        CASE('C     ')
!
        DO i=1,npoin
          x(i) = c
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('0     ')
!
        DO i=1,npoin
          x(i) = 0.d0
          IF (PRESENT (x_err)) THEN
            x_err(i)=0.d0
          ENDIF
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('Y     ')
!
        DO i=1,npoin
          x(i) = y(i)
          IF (PRESENT (x_err) .AND. PRESENT (y_err))THEN
            x_err(i) = y_err(i)
          ENDIF
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('+Y    ')
!
        DO i=1,npoin
          x(i) = y(i)
          IF (PRESENT (x_err).AND. PRESENT (y_err)) THEN
            x_err(i)=y_err(i)
          ENDIF
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('-Y    ')
!
        DO i=1,npoin
          x(i) = - y(i)
          IF (PRESENT (x_err).AND. PRESENT (y_err)) THEN
            x_err(i)=-y_err(i)
          ENDIF
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('1/Y   ')
!
        DO i=1,npoin
          x(i) = 1.d0/y(i)
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('Y+Z   ')
!
        DO i=1,npoin
          CALL twosum(y(i),z(i),x(i),error)
          IF (PRESENT (x_err) .AND. PRESENT (y_err)
     &      .AND. PRESENT (z_err)) THEN
            x_err(i)=(y_err(i)+z_err(i))+error
          ENDIF
          !X(I) = Y(I) + Z(I)
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('Y-Z   ')
!
        DO i=1,npoin
          CALL twosum(y(i),-z(i),x(i),error)
          IF (PRESENT (x_err) .AND. PRESENT (y_err)
     &      .AND. PRESENT (z_err)) THEN
            x_err(i)=(y_err(i)-z_err(i))+error
          ENDIF
          !X(I) = Y(I) - Z(I)
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('YZ    ')
!
        DO i=1,npoin
          CALL twoprod(y(i),z(i),x(i),error)
          IF (PRESENT (x_err) .AND. PRESENT (y_err)
     &      .AND. PRESENT (z_err)) THEN
            x_err(i)=(y(i) * z_err(i))+(y_err(i) * z(i))
     &               +(y_err(i) * z_err(i))
            x_err(i)=x_err(i)+ error
          ENDIF
          !X(I) = Y(I) * Z(I)
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('-YZ   ')
!
        DO i=1,npoin
          CALL twoprod(-y(i),z(i),x(i),error)
          IF (PRESENT (x_err) .AND. PRESENT (y_err)
     &      .AND. PRESENT (z_err)) THEN
            x_err(i)=(-y(i) * z_err(i))+(-y_err(i) * z(i))
            x_err(i)=x_err(i)+ error
          ENDIF
          !X(I) = -Y(I) * Z(I)
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('XY    ')
!
        DO i=1,npoin
          test=0.d0
          CALL twoprod(x(i),y(i),test,error)
          IF (PRESENT (x_err) .AND. PRESENT (y_err)) THEN
            x_err(i)=(x(i) * y_err(i))+(x_err(i) * y(i))
     &               +(x_err(i) * y_err(i))
            x_err(i)=x_err(i)+ error
          ENDIF
          x(i) = test
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('X+YZ  ')
!
        DO i=1,npoin
          CALL twoprod(y(i),z(i),temp(i),error1)
          tmp = x(i)
          CALL twosum(tmp,temp(i),x(i),error)
          IF (PRESENT (x_err).AND. PRESENT (y_err)
     &       .AND. PRESENT (z_err))THEN
            x_err(i)=x_err(i)+(y(i) * z_err(i))
     &               +(y_err(i) * z(i))
            x_err(i)=x_err(i)+ error1
            x_err(i)=x_err(i)+error
          ENDIF
          !X(I) = X(I) + Y(I) * Z(I)
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('X-YZ  ')
!
        DO i=1,npoin
          CALL twoprod(y(i),z(i),temp(i),error1)
          tmp = x(i)
          CALL twosum(tmp,-temp(i),x(i),error)
          IF (PRESENT (x_err).AND. PRESENT (y_err)
     &       .AND. PRESENT (z_err))THEN
            temp(i)=(y(i) * z_err(i))
     &              +(y_err(i) * z(i))+(y_err(i) * z_err(i))
            temp(i)=temp(i)+error1
            x_err(i)=(x_err(i)-temp(i))+error
          ENDIF
          !X(I) = X(I) - Y(I) * Z(I)
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('CXY   ')
!
        DO i=1,npoin
          x(i) = c * x(i) * y(i)
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('CYZ   ')
!
        DO i=1,npoin
          CALL twoprod(c*y(i),z(i),x(i),error)
          IF (PRESENT (x_err) .AND. PRESENT (y_err)
     &      .AND. PRESENT (z_err)) THEN
            x_err(i)= (c*y(i) * z_err(i))
     &               +(c*y_err(i) * z(i))+(c*y_err(i) * z_err(i))
            x_err(i)=  x_err(i)+error
          ENDIF
          !X(I) = C* Y(I) * Z(I)
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('CXYZ  ')
!
        DO i=1,npoin
          x(i) = c * x(i) * y(i) * z(i)
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('X+CYZ ')
!
        DO i=1,npoin
          x(i) = x(i) + c * y(i) * z(i)
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('Y/Z   ')
!
        DO i=1,npoin
          x(i) = y(i) / z(i)
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('CY/Z  ')
!
        DO i=1,npoin
          x(i) = c*y(i) / z(i)
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('CXY/Z ')
!
        DO i=1,npoin
          x(i) = c*x(i)*y(i) / z(i)
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('X+CY/Z')
!
        DO i=1,npoin
          x(i) = x(i) + c * y(i) / z(i)
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('X+Y   ')
!
        DO i=1,npoin
          tmp = x(i)
          CALL twosum(tmp,y(i),x(i),error)
          IF (PRESENT (x_err) .AND. PRESENT (y_err))THEN
            x_err(i)=(x_err(i)+y_err(i))+error
          ENDIF
          !X(I) = X(I) + Y(I)
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('X-Y   ')
        DO i=1,npoin
          tmp = x(i)
          CALL twosum(tmp,-y(i),x(i),error)
          IF (PRESENT (x_err))THEN
            x_err(i)=(x_err(i)-y_err(i))+error
          ENDIF
          !X(I) = X(I) - Y(I)
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('CX    ')
!
        DO i=1,npoin
          x(i) = c * x(i)
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('CY    ')
!
        DO i=1,npoin
          CALL twoprod(c,y(i),x(i),error)
          IF (PRESENT (x_err) .AND. PRESENT (y_err)) THEN
            x_err(i)=(c * y_err(i))
            x_err(i)=x_err(i)+error
          ENDIF
          !X(I) = C * Y(I)
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('Y+CZ  ')
!
        DO i=1,npoin
          x(i) = y(i) + c * z(i)
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('X+CY  ')
!
        DO i=1,npoin
          CALL twoprod(c,y(i),temp(i),error1)
          tmp = x(i)
          CALL twosum(tmp,temp(i),x(i),error)
          IF (PRESENT (x_err).AND. PRESENT (y_err))THEN
            x_err(i)=x_err(i)+(c* y_err(i))
            error=error + error1
            x_err(i)=x_err(i)+error
          ENDIF
          !X(I) = X(I) + C * Y(I)
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('SQR(Y)')
!
        DO i=1,npoin
          x(i) = sqrt(y(i))
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('ABS(Y)')
!
        DO i=1,npoin
          x(i) = abs(y(i))
          IF (PRESENT (x_err) .AND. PRESENT (y_err))THEN
            x_err(i) = abs(y_err(i))
          ENDIF
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('N(Y,Z)')
!
        DO i=1,npoin
          x(i) = sqrt( y(i)**2 + z(i)**2 )
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('Y+C   ')
!
        DO i=1,npoin
          x(i) = y(i) + c
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('X+C   ')
!
        DO i=1,npoin
          x(i) = x(i) + c
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('Y**C  ')
!
        DO i=1,npoin
          IF(y(i).GE.0.d0) THEN
            x(i) = y(i)**c
          ELSE
            WRITE(lu,101)
101         FORMAT(1x,'OV (BIEF): Y**C FORBIDDEN IF Y < 0')
            CALL plante(1)
            stop
          ENDIF
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('COS(Y)')
!
        DO i=1,npoin
          x(i) = cos(y(i))
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('SIN(Y)')
!
        DO i=1,npoin
          x(i) = sin(y(i))
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('ATN(Y)')
!
        DO i=1,npoin
          x(i) = atan(y(i))
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('A(Y,Z)')
!
        DO i=1,npoin
          x(i) = atan2(y(i),z(i))
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('+(Y,C)')
!
        DO i=1,npoin
          x(i) = max(y(i),c)
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('-(Y,C)')
!
        DO i=1,npoin
          x(i) = min(y(i),c)
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('+(Y,Z)')
!
        DO i=1,npoin
          x(i) = max(y(i),z(i))
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('-(Y,Z)')
!
        DO i=1,npoin
          x(i) = min(y(i),z(i))
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('YIFZ<C')
!
        DO i=1,npoin
          IF ( z(i).LT.c ) x(i) = y(i)
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE('C(Y-Z)')
!
        DO i=1,npoin
          x(i) = c*(y(i)-z(i))
        ENDDO
!
!-----------------------------------------------------------------------
!
        CASE DEFAULT
!
          WRITE(lu,1001) op
1001      FORMAT(1x,'OV (BIEF) : UNKNOWN OPERATION: ',a8)
          CALL plante(1)
          stop
!
!-----------------------------------------------------------------------
!
      END SELECT
!
!-----------------------------------------------------------------------
!
      RETURN
      END