OperatorJ.f90

Go to the documentation of this file.

!*********************************************************************
!* Illinois Open Source License                                      *
!*                                                                   *
!* University of Illinois/NCSA                                       * 
!* Open Source License                                               *
!*                                                                   *
!* Copyright@2008, University of Illinois.  All rights reserved.     *
!*                                                                   *
!*  Developed by:                                                    *
!*                                                                   *
!*     Center for Simulation of Advanced Rockets                     *
!*                                                                   *
!*     University of Illinois                                        *
!*                                                                   *
!*     www.csar.uiuc.edu                                             *
!*                                                                   *
!* Permission is hereby granted, free of charge, to any person       *
!* obtaining a copy of this software and associated documentation    *
!* files (the "Software"), to deal with the Software without         *
!* restriction, including without limitation the rights to use,      *
!* copy, modify, merge, publish, distribute, sublicense, and/or      *
!* sell copies of the Software, and to permit persons to whom the    *
!* Software is furnished to do so, subject to the following          *
!* conditions:                                                       *
!*                                                                   *
!*                                                                   *
!* @ Redistributions of source code must retain the above copyright  * 
!*   notice, this list of conditions and the following disclaimers.  *
!*                                                                   * 
!* @ Redistributions in binary form must reproduce the above         *
!*   copyright notice, this list of conditions and the following     *
!*   disclaimers in the documentation and/or other materials         *
!*   provided with the distribution.                                 *
!*                                                                   *
!* @ Neither the names of the Center for Simulation of Advanced      *
!*   Rockets, the University of Illinois, nor the names of its       *
!*   contributors may be used to endorse or promote products derived * 
!*   from this Software without specific prior written permission.   *
!*                                                                   *
!* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,   *
!* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   *
!* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          *
!* NONINFRINGEMENT.  IN NO EVENT SHALL THE CONTRIBUTORS OR           *
!* COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER       * 
!* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,   *
!* ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE    *
!* USE OR OTHER DEALINGS WITH THE SOFTWARE.                          *
!*********************************************************************
!* Please acknowledge The University of Illinois Center for          *
!* Simulation of Advanced Rockets in works and publications          *
!* resulting from this software or its derivatives.                  *
!*********************************************************************
SUBROUTINE operatorj (Lm, Ld, Lb, L2, Lbar,Mm, M2, cm, cb, cd, &
     dbm, dbb, dbd, dmm,dmb,dmd, ddm, ddb, ddd, &
     am, ab, ad, a2, km, k2, gm, g2, vm, v2, a_eta, a_zeta, &  
     opj, n, ninv, r, rinv, q2, qm, w)

  IMPLICIT NONE

  integer :: i

  REAL*8 :: cm, cb, cd
  REAL*8 :: a_eta, a_zeta

  REAL*8 :: gm, g2, km, k2, vm, v2

  REAL*8, DIMENSION(1:6,1:6) :: am, ab, ad, a2
  REAL*8, DIMENSION(1:6,1:6) :: lm, ld, lb, l2, lbar
  REAL*8, DIMENSION(1:6,1:6) :: mm, m2
  REAL*8, DIMENSION(1:6,1:6) :: dbm, dbb, dbd, dmm,dmb,dmd, ddm, ddb, ddd
!                                                       *              *
  
  REAL*8, DIMENSION(1:6,1:6)  :: opj , test
  
  REAL*8, DIMENSION(1:6,1:6) :: n, ninv, r, w, rinv

  REAL*8, DIMENSION(1:6,1:6) :: atmp, lr_lbarinv
  REAL*8, DIMENSION(1:6,1:6) :: bm, bb, bd
  REAL*8, DIMENSION(1:6,1:6) :: qm, qb, qd, q2

  REAL*8, DIMENSION(1:6,1:6) :: ident = reshape( &
       (/1.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, &
         0.d0, 1.d0, 0.d0, 0.d0, 0.d0, 0.d0, &
         0.d0, 0.d0, 1.d0, 0.d0, 0.d0, 0.d0, &
         0.d0, 0.d0, 0.d0, 1.d0, 0.d0, 0.d0, &
         0.d0, 0.d0, 0.d0, 0.d0, 1.d0, 0.d0, &
         0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 1.d0 /),(/6,6/) )

! W


  CALL constant_w(gm, g2, k2, km, a_eta, a_zeta, vm, v2, w)

!  -1
! N  

  n = matmul(l2,ddd - ident) - matmul(w,matmul(lm,dmd))

  CALL invert2(n, ninv, 6)

! R

  r = matmul(w,matmul(lm,am)) - matmul(l2,a2)

  CALL invert2(r,rinv,6)

! -- Q --

! Qm
  atmp = lm - lbar
  CALL invert2( atmp, lr_lbarinv, 6)

  qm = matmul((ident - am),lr_lbarinv)

!Q2

  atmp = l2 - lbar
  CALL invert2( atmp, lr_lbarinv, 6)
  
  q2 = matmul(ident - a2,lr_lbarinv)

! Qb

  qb = q2

! Qd

  qd = q2

!                   -1  
! Bm = [ ( Dmd - I)N  (L2 Q2-W Lm Qm )-Qm } A2L2NR

!test  = MATMUL( Dmd * Ninv*(L2*Q2 - W*Lm*Qm) - Qm,TRANSPOSE(A2)*L2*Ninv*R)
!!$  Bb = ( Dbd * Ninv*(L2*Q2 - W*Lm*Qm) - Qb)*TRANSPOSE(A2)*L2*Ninv*R
!!$  Bd = ( (Ddd - Ident) * Ninv*(L2*Q2 - W*Lm*Qm) - Qd)*TRANSPOSE(A2)*L2*Ninv*R



!!$  Bm = MATMUL( MATMUL (Dmd , MATMUL(Ninv,( MATMUL(L2,Q2) - MATMUL(W,MATMUL(Lm,Qm)) ) ) ) - Qm, &
!!$       MATMUL(TRANSPOSE(A2),MATMUL(L2,MATMUL(Ninv,R))) )


  bm = matmul(l2,q2) - matmul(w, matmul(lm,qm) )

  bm = matmul(ninv,bm)

  bm = matmul(dmd,bm)

  bm = bm - qm

  bm = matmul(bm,transpose(a2))
  
  bm = matmul(bm,l2)

  bm = matmul(bm,ninv)

  bm = matmul(bm,r)

!!$  DO i = 1,6 
!!$     PRINT*,'now',Bm(i,1:6)
!!$     PRINT*,'org',test(i,1:6)
!!$  ENDDO
  
  bb = matmul(l2,q2) - matmul(w, matmul(lm,qm) )

  bb = matmul(ninv,bb)

  bb = matmul(dbd,bb)

  bb = bb - qb

  bb = matmul(bb,transpose(a2))
  
  bb = matmul(bb,l2)

  bb = matmul(bb,ninv)

  bb = matmul(bb,r)

!!$  Bb = MATMUL( MATMUL( Dbd , MATMUL(Ninv,( MATMUL(L2,Q2) - MATMUL(W,MATMUL(Lm,Qm)) ) ) ) - Qb , &
!!$       MATMUL(TRANSPOSE(A2),MATMUL(L2,MATMUL(Ninv,R))) )


  bd = matmul(l2,q2) - matmul(w, matmul(lm,qm) )

  bd = matmul(ninv,bd)

  bd = matmul(ddd - ident,bd)

  bd = bd - qd

  bd = matmul(bd,transpose(a2))
  
  bd = matmul(bd,l2)

  bd = matmul(bd,ninv)

  bd = matmul(bd,r)

!!$  Bd = MATMUL( MATMUL((Ddd - Ident), MATMUL(Ninv, (MATMUL(L2,Q2) - MATMUL(W,MATMUL(Lm,Qm))) ) ) - Qd, &
!!$       MATMUL(TRANSPOSE(A2),MATMUL(L2,MATMUL(Ninv,R))) )

! OpJ = cm*Lm*Bm + cb*L2*Bb + Cd*L2*Bd + L2*( (Ddd - Ident) - Dbd)* Ninv * R

  opj = matmul( ( ddd - ident) - dbd, ninv)

  opj = matmul( l2, opj)

  opj = matmul( opj, r)

  
  opj = cm*matmul(lm,bm) + cb*matmul(l2,bb) + cd*matmul(l2,bd) + opj



!!$  OpJ = cm*MATMUL(Lm,Bm) + cb*MATMUL(L2,Bb) + cd*MATMUL(L2,Bd) + &
!!$       MATMUL(L2, MATMUL( ( Ddd - Ident) - Dbd, MATMUL(Ninv , R)) )

!!$  DO i = 1, 6
!!$     PRINT*,OpJ(i,1:6)
!!$  ENDDO

END SUBROUTINE operatorj