C  PCPRB2.FOR  27 February 1991
C
      PROGRAM PCPRB2
C
C  The aircraft stability problem.
C
C  The variables X(I) are control parameters for an aircraft under
C  the guidance of a pilot:
C
C
C  X(1) = Roll
C  X(2) = Pitch
C  X(3) = Yaw
C  X(4) = Angle of attack
C  X(5) = Sideslip
C  X(6) = Elevator
C  X(7) = Aileron
C  X(8) = Rudder
C
C
C  The function F is of the following form:
C
C
C  For I=1 to 5,
C
C  F(I)=Sum(J=1,8) Sum(K=1,8) B(I,J)*X(J)+PHI(I,J,K)*X(J)*X(K)
C
C  F(6)=X(IFIX1)-VAL1
C  F(7)=X(IFIX2)-VAL2
C
C
C  Options:
C
C
C  ICHOOZ=1  IFIX1=6, VAL1=-.050, Rheinboldt BARRAY.
C  ICHOOZ=2  IFIX1=6, VAL1=-.008, Rheinboldt BARRAY.
C  ICHOOZ=3  IFIX1=6, VAL1=0.000, Rheinboldt BARRAY.
C  ICHOOZ=4  IFIX1=6, VAL1=0.050, Rheinboldt BARRAY.
C  ICHOOZ=5  IFIX1=6, VAL1=0.100, Rheinboldt BARRAY.
C  ICHOOZ=6  IFIX1=6, VAL1=-0.01250, Rheinboldt BARRAY.  (Bifurcation)
C
C  ICHOOZ=11  IFIX1=6, VAL1=-.050, Melhem BARRAY.
C  ICHOOZ=12  IFIX1=6, VAL1=-.008, Melhem BARRAY.
C  ICHOOZ=13  IFIX1=6, VAL1=0.000, Melhem BARRAY.
C  ICHOOZ=14  IFIX1=6, VAL1=0.050, Melhem BARRAY.
C  ICHOOZ=15  IFIX1=6, VAL1=0.100, Melhem BARRAY.
C  ICHOOZ=16  IFIX1=6, VAL1=-0.01250, Melhem BARRAY.  (Bifurcation)
C
C  For all choices of ICHOOZ, IFIX2=8, VAL2=0.0.
C  If ICHOOZ=6 or 16, LIM=0, otherwise LIM=7.
C
C
C  Limit points:
C
C
C  Melhem lists the following limit points in X7.  Note that Melhem
C  has BARRAY(4,1)=1.0, BARRAY(4,2)=0.0:
C
C
C     X1       X2        X3        X4        X5      X6     X7      X8
C
C  -2.9691  0.83074  -0.072748  0.41029  -0.26880  -.05  0.50919    0.0
C  -2.8158 -0.17481  -0.089469  0.026319  0.070951 -.008 0.20442    0.0
C  -3.7571 -0.64911  -0.39350   0.091810  0.19685  -.008 -.0038238  0.0
C  -4.1637  0.092284 -0.092610  0.022402 -0.017106 -.008 0.37823    0.0
C  -2.5839 -0.22128  -0.054079  0.013524  0.090871 0.0   0.18608    0.0
C  -3.9007 -1.1421   -0.57863   0.13284   0.32685  0.0  -0.50703    0.0
C  -2.3610 -0.72360   0.032739 -0.039108  0.29347  0.05  0.29272    0.0
C  -2.2982  1.4033    0.063244 -0.079383  0.58336  0.10  0.58336    0.0
C
C
C  Rheinboldt lists the following limit points, using
C  BARRAY(4,1)=0.0, BARRAY(4,2)=1.0:
C
C
C     X1       X2        X3        X4        X5      X6     X7      X8
C
C   2.9648  0.82556   0.07366   0.041309  0.26734 -0.050 -.050481   0.0
C   2.8173 -0.17628   0.08992   0.026429 -0.07147 -0.008 -.204973   0.0
C   3.7579 -0.65541   0.38658   0.092520 -0.19867 -0.008 0.006200   0.0
C   4.1638  0.08913   0.09480   0.022888  0.16232 -0.008 -.377660   0.0
C   2.5873 -0.22354   0.05468   0.013676 -0.09168  0.000 -.186908   0.0
C   3.9005 -1.14815   0.58156   0.133516 -0.32858   0.000  .510158  0.0
C   2.3639 -0.72974  -0.31604  -0.038785 -0.29583   0.050 -.295772  0.0
C   2.2992 -1.41023  -0.06184  -0.079009 -0.58629   0.100 -.689717  0.0
C
C
C  Bifurcation points:
C
C
C  Rheinboldt, using BARRAY(4,1)=0.0, BARRAY(4,2)=1.0, lists:
C
C
C   4.482   0.1632    0.02373   0.006205  0.03527 -0.0006177 -0.3986 0.0
C   3.319  -0.1869    0.1605    0.04379  -0.06888 -0.01250   -0.2374 0.0
C   4.466   0.1467    0.04045   0.009777  0.03089 -0.006129  -0.3995 0.0
C  -3.325   0.1880   -0.1614    0.04395   0.06911 -0.01247    0.2367 0.0
C
C  References:
C
C  A A Schy, M E Hannah,
C  Prediction of Jump Pheonomena in Roll-Coupled Maneuvers of Airplanes,
C  Journal of Aircraft, 14, 1977,  Pages 375-382.
C
C  J W Young, A A Schy, K G Johnson
C  Prediction of Jump Phenomena in Aircraft Maneuvers, Including Nonlinear
C    Aerodynamic Effects,
C  Journal of Guidance and Control 1, 1978, Pages 26-31.
C
      INTEGER   LIW
      INTEGER   LRW
      INTEGER   NVAR
C
      PARAMETER (NVAR=8)
      PARAMETER (LIW=29+NVAR)
      PARAMETER (LRW=29+NVAR*(NVAR+4))
C
      EXTERNAL  DENSLV
      EXTERNAL  FPAIR
      EXTERNAL  FXAIR
      EXTERNAL  PITCON
C
      REAL      BARRAY(5,8)
      REAL      HMAX
      INTEGER   I
      INTEGER   ICHOOZ
      INTEGER   IERROR
      INTEGER   IFIX1
      INTEGER   IFIX2
      INTEGER   IPAR(2)
      INTEGER   IWORK(LIW)
      INTEGER   J
      INTEGER   K
      INTEGER   LIM
      INTEGER   LOUNIT
      CHARACTER NAME*12
      REAL      RPAR(42)
      REAL      RWORK(LRW)
      REAL      VAL1
      REAL      VAL2
      REAL      XR(NVAR)
C
      LOUNIT=6
C
C  Set options.
C
      ICHOOZ=1
      IFIX1=6
      IFIX2=8
C
C  Initialize work arrays.
C
      DO 4 I=1,LIW
        IWORK(I)=0
4       CONTINUE
      DO 5 I=1,LRW
        RWORK(I)=0.0
5       CONTINUE
C
C  Set VAL1, VAL2 based on ICHOOZ
C
      IF(ICHOOZ.EQ.1.OR.ICHOOZ.EQ.11)VAL1=-.050
      IF(ICHOOZ.EQ.2.OR.ICHOOZ.EQ.12)VAL1=-.008
      IF(ICHOOZ.EQ.3.OR.ICHOOZ.EQ.13)VAL1=0.000
      IF(ICHOOZ.EQ.4.OR.ICHOOZ.EQ.14)VAL1=0.050
      IF(ICHOOZ.EQ.5.OR.ICHOOZ.EQ.15)VAL1=0.100
      IF(ICHOOZ.EQ.6.OR.ICHOOZ.EQ.16)VAL1=-0.01250
      VAL2=0.0
C
C  Set some parameters
C
      IPAR(1)=IFIX1
      IPAR(2)=IFIX2
      IF(ICHOOZ.NE.6.AND.ICHOOZ.NE.16)THEN
        HMAX=0.4
      ELSE
        HMAX=0.2
        ENDIF
      IF(ICHOOZ.NE.6.AND.ICHOOZ.NE.16)THEN
        LIM=7
      ELSE
        LIM=0
        ENDIF
C
C  Set starting point
C
      DO 10 I=1,NVAR
        XR(I)=0.0
10      CONTINUE
      XR(IFIX1)=VAL1
      XR(IFIX2)=VAL2
C
C  Set BARRAY
C
      BARRAY(1,1)=-3.933
      BARRAY(2,1)=0.0
      BARRAY(3,1)=0.002
      IF(ICHOOZ.LT.10)THEN
        BARRAY(4,1)=0.0
      ELSE
        BARRAY(4,1)=1.0
        ENDIF
      BARRAY(5,1)=0.0
      BARRAY(1,2)=0.107
      BARRAY(2,2)=-0.987
      BARRAY(3,2)=0.0
      IF(ICHOOZ.LT.10)THEN
        BARRAY(4,2)=1.0
      ELSE
        BARRAY(4,2)=0.0
        ENDIF
      BARRAY(5,2)=0.0
      BARRAY(1,3)=0.126
      BARRAY(2,3)=0.0
      BARRAY(3,3)=-0.235
      BARRAY(4,3)=0.0
      BARRAY(5,3)=-1.0
      BARRAY(1,4)=0.0
      BARRAY(2,4)=-22.95
      BARRAY(3,4)=0.0
      BARRAY(4,4)=-1.0
      BARRAY(5,4)=0.0
      BARRAY(1,5)=-9.99
      BARRAY(2,5)=0.0
      BARRAY(3,5)=5.67
      BARRAY(4,5)=0.0
      BARRAY(5,5)=-0.196
      BARRAY(1,6)=0.0
      BARRAY(2,6)=-28.37
      BARRAY(3,6)=0.0
      BARRAY(4,6)=-0.168
      BARRAY(5,6)=0.0
      BARRAY(1,7)=-45.83
      BARRAY(2,7)=0.0
      BARRAY(3,7)=-0.921
      BARRAY(4,7)=0.0
      BARRAY(5,7)=-0.0071
      BARRAY(1,8)=-7.64
      BARRAY(2,8)=0.0
      BARRAY(3,8)=-6.51
      BARRAY(4,8)=0.0
      BARRAY(5,8)=0.0
C
C  Copy BARRAY, VAL1, VAL2 into RPAR.
C
      K=0
      DO 30 I=1,5
        DO 29 J=1,NVAR
          K=K+1
          RPAR(K)=BARRAY(I,J)
29        CONTINUE
30      CONTINUE
      K=K+1
      RPAR(K)=VAL1
      K=K+1
      RPAR(K)=VAL2
C
C  Set work array input
C
C  IWORK(1)=0 ; This is a startup
C  IWORK(2)=7 ; Use X(7) for initial parameter
C  IWORK(3)=0 ; Program may choose parameter index
C  IWORK(4)=0 ; Update jacobian every newton step
C  IWORK(5)=0 ; Seek no target values.
C  IWORK(6)=LIM ; Seek limit points in X(LIM)
C  IWORK(7)=1 ; Control amount of output.
C  IWORK(8)=6 ; Output unit
C  IWORK(9)=0 ; Jacobian choice.
C
      IWORK(1)=0
      IWORK(2)=7
      IWORK(3)=0
      IWORK(4)=0
      IWORK(5)=0
      IWORK(6)=LIM
      IWORK(7)=1
      IWORK(8)=LOUNIT
      IWORK(9)=0
C
C  RWORK(1)=0.0001 ; Absolute error tolerance
C  RWORK(2)=0.0001 ; Relative error tolerance
C  RWORK(3)=0.0001 ; Minimum stepsize
C  RWORK(4)=HMAX   ; Maximum stepsize
C  RWORK(5)=0.1    ; Starting stepsize
C  RWORK(6)=-1.0   ; Starting direction
C  RWORK(7)=0.0    ; Target value
C
      RWORK(1)=0.0001
      RWORK(2)=0.0001
      RWORK(3)=0.0001
      RWORK(4)=HMAX
      RWORK(5)=0.1
      RWORK(6)=-1.0
      RWORK(7)=0.0
C
      WRITE(6,*)' '
      WRITE(6,*)'PCPRB2:'
      WRITE(6,*)'PITCON sample program.'
      WRITE(6,*)' '
      WRITE(6,*)'The Aircraft Stability Problem'
      WRITE(6,*)'Using option ',ICHOOZ
      WRITE(6,*)'Fix variable ',IFIX1,' at ',VAL1
      WRITE(6,*)'Fix variable ',IFIX2,' at ',VAL2
      IF(ICHOOZ.LT.10)THEN
        WRITE(6,*)'Rheinboldt version of BARRAY'
      ELSE
        WRITE(6,*)'Melhem version of BARRAY'
        ENDIF
      WRITE(6,*)' '
      WRITE(6,*)'Step  Type of Point     '
      WRITE(6,*)' '
      I=0
      NAME='Start point  '
      WRITE(6,'(1X,I3,2X,A12,2X,8F7.3)')I,NAME,(XR(J),J=1,NVAR)
C
      DO 40 I=1,50
        CALL PITCON(FPAIR,RPAR,FXAIR,IERROR,IPAR,IWORK,LIW,
     *  NVAR,RWORK,LRW,XR,DENSLV)
C
        IF(IWORK(1).EQ.1)THEN
          NAME='Corrected    '
        ELSEIF(IWORK(1).EQ.2)THEN
          NAME='Continuation '
        ELSEIF(IWORK(1).EQ.3)THEN
          NAME='Target point '
        ELSEIF(IWORK(1).EQ.4)THEN
          NAME='Limit point  '
          ENDIF
C
        WRITE(6,'(1X,I3,2X,A12,2X,8F7.3)')I,NAME,(XR(J),J=1,NVAR)
C
        IF(IERROR.NE.0)THEN
          WRITE(6,*)'Iteration halted, IERROR=',IERROR
          GO TO 50
          ENDIF
C
40      CONTINUE
50    CONTINUE
C
      WRITE(6,*)' '
      WRITE(6,*)'Jacobians:     ',IWORK(19)
      WRITE(6,*)'Factorizations:',IWORK(20)
      WRITE(6,*)'Solves:        ',IWORK(21)
      WRITE(6,*)'Functions:     ',IWORK(22)
C
      STOP
      END
      SUBROUTINE FXAIR(NVAR,RPAR,IPAR,X,FX,IERROR)
C
C***********************************************************************
C
      INTEGER   NVAR
C
      REAL      FX(NVAR)
      INTEGER   I
      INTEGER   IERROR
      INTEGER   IFIX1
      INTEGER   IFIX2
      INTEGER   IPAR(2)
      INTEGER   J
      INTEGER   K
      REAL      PHI
      REAL      RPAR(42)
      REAL      VAL1
      REAL      VAL2
      REAL      X(NVAR)
C
C  Compute linear terms
C
      K=0
      DO 20 I=1,5
        FX(I)=0.0
        DO 10 J=1,8
          K=K+1
          FX(I)=FX(I)+RPAR(K)*X(J)
10        CONTINUE
20      CONTINUE
C
C  Compute nonlinear terms
C
      PHI=-.727*X(2)*X(3)+8.39*X(3)*X(4)-684.4*X(4)*X(5)+63.5*X(4)*X(7)
      FX(1)=FX(1)+PHI
C
      PHI=.949*X(1)*X(3)+.173*X(1)*X(5)
      FX(2)=FX(2)+PHI
C
      PHI=-.716*X(1)*X(2)-1.578*X(1)*X(4)+1.132*X(4)*X(7)
      FX(3)=FX(3)+PHI
C
      PHI=-X(1)*X(5)
      FX(4)=FX(4)+PHI
C
      PHI=X(1)*X(4)
      FX(5)=FX(5)+PHI
C
C  Two function values restrict two variables:
C
      IFIX1=IPAR(1)
      VAL1=RPAR(41)
      FX(6)=X(IFIX1)-VAL1
C
      IFIX2=IPAR(2)
      VAL2=RPAR(42)
      FX(7)=X(IFIX2)-VAL2
C
      RETURN
      END
      SUBROUTINE FPAIR(NVAR,RPAR,IPAR,X,FPRIME,IERROR)
C
C***********************************************************************
C
      INTEGER   NVAR
C
      REAL      FPRIME(NVAR,NVAR)
      INTEGER   I
      INTEGER   IERROR
      INTEGER   IFIX1
      INTEGER   IFIX2
      INTEGER   IPAR(2)
      INTEGER   J
      INTEGER   K
      REAL      RPAR(42)
      REAL      X(NVAR)
C
C  Linear part of F
C
      K=0
      DO 20 I=1,5
        DO 10 J=1,8
          K=K+1
          FPRIME(I,J)=RPAR(K)
10        CONTINUE
20      CONTINUE
C
      DO 40 I=6,7
        DO 30 J=1,8
          FPRIME(I,J)=0.0
30        CONTINUE
40      CONTINUE
C
      IFIX1=IPAR(1)
      IFIX2=IPAR(2)
      FPRIME(6,IFIX1)=1.0
      FPRIME(7,IFIX2)=1.0
C
C  Nonlinear part of F
C
      FPRIME(1,2)=FPRIME(1,2)-.727*X(3)
      FPRIME(1,3)=FPRIME(1,3)-.727*X(2)+8.39*X(4)
      FPRIME(1,4)=FPRIME(1,4)+8.39*X(3)-684.4*X(5)+63.5*X(7)
      FPRIME(1,5)=FPRIME(1,5)-684.4*X(4)
      FPRIME(1,7)=FPRIME(1,7)+63.5*X(4)
C
      FPRIME(2,1)=FPRIME(2,1)+.949*X(3)+.173*X(5)
      FPRIME(2,3)=FPRIME(2,3)+.949*X(1)
      FPRIME(2,5)=FPRIME(2,5)+.173*X(1)
C
      FPRIME(3,1)=FPRIME(3,1)-.716*X(2)-1.578*X(4)
      FPRIME(3,2)=FPRIME(3,2)-.716*X(1)
      FPRIME(3,4)=FPRIME(3,4)-1.578*X(1)+1.132*X(7)
      FPRIME(3,7)=FPRIME(3,7)+1.132*X(4)
C
      FPRIME(4,1)=FPRIME(4,1)-X(5)
      FPRIME(4,5)=FPRIME(4,5)-X(1)
C
      FPRIME(5,1)=FPRIME(5,1)+X(4)
      FPRIME(5,4)=FPRIME(5,4)+X(1)
C
      RETURN
      END