$debug
      use msflib
      logical(4) status
      external rawfile, pltfile,go,quad,sig,sig2,prst,least,
     >         leasthelp,regre,regrehelp,inst,insthelp
      status=aboutboxqq
     >('ev Copyright Y. Fujii, Jan. 15, 1998'c)
*
      status=insertmenuqq(3,0,$menuenabled,'Instron'c,nul)
*      status=appendmenuqq(3,$menuenabled,'Raw file'c,rawfile)
*      status=appendmenuqq(3,$menuenabled,'Plt file'c,pltfile)
      status=appendmenuqq(3,$menuenabled,'Go!'c,inst)
      status=appendmenuqq(3,$menuenabled,'help'c,insthelp)
*
      status=insertmenuqq(4,0,$menuenabled,'Sig 2-D'c,nul)
      status=appendmenuqq
     >(4,$menuenabled,'σx,σy,τxy→σ1,σ2,α'c,sig)
      status=appendmenuqq
     >(4,$menuenabled,
     >'σx,σy,τxy,α→σx'',σy'',τx''y'''c,sig2)
*
      status=insertmenuqq(5,0,$menuenabled,'Sig 3-D'c,nul)
      status=appendmenuqq(5,$menuenabled,
     >'σx,σy,σz,τxy,τyz,τzx→σ1,σ2,σ3,φ,θ'c,prst)
*
      status=insertmenuqq(6,0,$menuenabled,'Quadratic equation'c,nul)
      status=appendmenuqq(6,$menuenabled,
     >'ax2+bx+c=0→x1, x2'c,quad)
*
      status=insertmenuqq(7,0,$menuenabled,'Least squares'c,nul)
      status=appendmenuqq(7,$menuenabled,'Go!'c,least)
*      status=appendmenuqq(7,$menuenabled,'Go!'c,least)
      status=appendmenuqq(7,$menuenabled,'Instruction'c,leasthelp)
*
      status=insertmenuqq(8,0,$menuenabled,'Regression line'c,nul)
*      status=appendmenuqq(8,$menuenabled,'Data file'c,rawfile)
      status=appendmenuqq(8,$menuenabled,'Go!'c,regre)
      status=appendmenuqq(8,$menuenabled,'Instruction'c,regrehelp)
*
      open(7,file='user')
      write(7,*)'Welcome to ev!'
      write(7,*)'Choice Instron, Sig 2-D, Sig 3-D, Quadratic equation,'
      write(7,*)'Least squares, or Regression line.'
      iresult=incharqq()
      end
*
      subroutine regrehelp(status)
      use msflib
      logical status
      open(7,file='user')
      write(7,*)'これは重みなし最小二乗法を用いて回帰直線の'
      write(7,*)'傾きとy切片の最確値と標準偏差を求めるプログラムです。'
      write(7,*)'通常は相関係数を評価しますが、'
      write(7,*)'x値に対してy値に変動がない場合にはデータが'
      write(7,*)'直線になっていても極めて小さい相関係数が評価されます。'
      write(7,*)'これは、ばらつきが少なくても「相関」がないからです。'
      write(7,*)'このような時にデータのばらつきの度合いを評価するために'
      write(7,*)'使って下さい。最確値は当然回帰直線と同じです。'
      write(7,*)'Data fileを指定してからGo!をクリックして下さい。'
      write(7,*)'今のところ最大500データに対応しています。'
      write(7,*)'データファイルにはx値,y値を1行に1組づつ書いて下さい。'
      return
      end
*
      subroutine insthelp(status)
      use msflib
      logical status
      open(7,file='user')
      write(7,*)'これはMerlinを使って計測したrawファイルから'
      write(7,*)'単軸圧縮強度、ヤング率、ポアソン比、限界ひずみを'
      write(7,*)'求め、時間、余計な変形を補正したひずみ、応力を'
      write(7,*)'出力するプログラムです。'
      return
      end
*progeam least
      subroutine regre(status)
      use msflib
      logical status
      DIMENSION A(500,2),B(500),ANS(2),ERR(2),
     >          AA(2,2),AB(2),AB2(2),AR(2,2)
      open(5,file=' ')
      open(7,file='user')
* MAX NLINE=500
      N1=500
      N2=2
*      open(5,file=' ')
      ncol=2
      nline=0
      DO 200 K=1,1000
      READ(5,*,end=9997)A(K,1),b(k)
      a(k,2)=1.
      nline=nline+1
  200 CONTINUE
 9997 continue
      close(5)
*  300 FORMAT(7F10.3)
*
      CALL LS(A,N1,N2,NLINE,NCOL,B,ANS,ERR,ILL1,ILL2,AA,AB,AB2,AR)
*
      write(7,*)'傾きの最確値と標準偏差'
      WRITE(7,1000)ANS(1),ERR(1)
      write(7,*)'y切片の最確値と標準偏差'
      WRITE(7,1000)ANS(2),ERR(2)
 1000 FORMAT(1H ,2G12.4)
  500 CONTINUE
      return
      END
*
      subroutine leasthelp(status)
      use msflib
      logical status
      open(7,file='user')
      write(7,*)'これは重みなし最小二乗法を用いて未知数の最確値と'
      write(7,*)'誤差（標準偏差）を求めるプログラムです。'
      write(7,*)'Data fileを指定してからGo!をクリックして下さい。'
      write(7,*)'今のところ最大7カラム、500行に対応しています。'
      write(7,*)'観測方程式をB=AXとして、'
      write(7,*)'データファイルのフォーマットは以下のようです。'
      write(7,*)'1行目：行数(nline)、カラム数(ncol)'
      write(7,*)'2行目：A(1,1),,, A(1,ncol), B(1)'
      write(7,*)'n行目：A(n-1,1),,,A(n-1,ncol), B(n-1)'
      write(7,*)'nline+1行目：A(nline,1),,,A(nline,ncol),B(nline)'
      return
      end
*progeam least
      subroutine least(status)
      use msflib
      logical status
      DIMENSION A(500,7),B(500),ANS(7),ERR(7),
     >          AA(7,7),AB(7),AB2(7),AR(7,7)
      open(7,file='user')
      open(5,file=' ')
*progeam least
*CODED BY Y.FUJII, DEC.3,1992
*
* MAX NCOL=7
* MAX NLINE=500
*
*  NLINE  NCOL (2I5)
*  A(1,1)------------A(1,NCOL)     B(1)
*    |                 |             |
*    |                 |             |
*  A(NLINE,1)--------A(NLINE,NCOL) B(NLINE)
*
      N1=500
      N2=7
*      open(5,file=' ')
      READ(5,*)NLINE,NCOL
*  100 FORMAT(2I5)
      DO 200 K=1,NLINE
      READ(5,*)(A(K,J),J=1,NCOL),b(k)
  200 CONTINUE
      close(5)
*  300 FORMAT(7F10.3)
*
      CALL LS(A,N1,N2,NLINE,NCOL,B,ANS,ERR,ILL1,ILL2,AA,AB,AB2,AR)
*
      write(7,*)'番号・最確値・標準偏差'
      DO 500 K=1,NCOL
      WRITE(7,1000)k,ANS(K),ERR(K)
 1000 FORMAT(1H ,I5, 2G12.4)
  500 CONTINUE
      return
      END

*
* solve quadratic eauation
*
      subroutine quad(status)
      use msflib
      logical status
*      write(*,*)'Solve ax^2+bx+c=0'
      open(7,file='user')
      write(7,*)'Input a, b, c'
      read(7,*)a,b,c
      det=b**2-4.*a*c
      if(det.lt.0.)then
      write(7,*)'can''t solve'
      else
      x1=(-b+sqrt(det))/(2.*a)
      x2=(-b-sqrt(det))/(2.*a)
      write(7,*)'x1, x2',x1,x2
      endif
*      write(7,*)'Hit any key'
*      iresult=incharqq()
*      close(7)
      return
      end
*
* rawfile name
*
      subroutine inst(status)
      use msflib
      dimension fl(5000),str1(5000),str2(5000),time(5000)
      DIMENSION A(4999,2),B(4999),ANS(2),ERR(2),
     >          AA(2,2),AB(2),AB2(2),AR(2,2),disp(5000),strrr(5000)
      logical status
      character*72 dumc
* ev(ver.2.1)
* create plot file
      open(5,file=' ')
      open(6,file=' ')

      N1=4999
      N2=2
      nnn=(n1-1)/2
*
      open(7,file='user')
      write(7,*)'供試体の直径と長さを「30,60」のようにmm単位で入力し'
      write(7,*)'て下さい。'
      read(7,*)dia,flen
*      close(7)
      area=dia**2*3.141593/4.
*
*      write(*,*)'input c1,c2'
*      read(*,*)c1,c2
      c1=0.4
      c2=0.6
      do 9100 ll=1,4
      read(5,2224)dumc
 2224 format(a72)
 9100 continue
*
*      r=14.75
*      area=r**2*3.141593
* size is in mm, load is in N
      do 100 k=1,10000
      read(5,*,end=9999)time(k),disp(k),str1(k),str2(k),fl(k)
  100 continue
 9999 continue
      flmax=10000.
* search maximum load
      do 200 j=1,k-1
      if(fl(j).le.flmax)then
      flmax=fl(j)
      kmax=j
      endif
  200 continue
      kend=k-1
* search loading start point
      do 233 j=1,k-1
      if(fl(j).le.-100.)then
      ksta=j-1+4
      dsta=-disp(j-1)
      goto234
      endif
  233 continue
*
  234 continue
      smax=-flmax/area
      ecc=-str1(kmax)
      etc=str2(kmax)
*      write(6,1000)smax,ecc,etc
* 1000 format(1h ,'smax(MPa),eac(%),elc(%)',f10.1,2f10.2)
*
      do 298 j40=1,k-1
      if(fl(j40).lt.flmax*c1)goto398
  298 continue
  398 continue
*
      do 300 j60=1,k-1
      if(fl(j60).lt.flmax*c2)goto400
  300 continue
  400 continue
 2222 continue
*      write(6,3000)j40,j60
 3000 format(1h ,'elastic constants are estimated data number',2i4)
*
      jjj=j60-j40+1
      nline=jjj
      ncol=2
*
      do 500 i=1,jjj
      k=i+j40-1
      a(i,1)=-str1(k)
      a(i,2)=1.
      b(i)=-fl(k)/area
*      write(6,*)a(i,1),a(i,2),b(i)
  500 continue
*
      CALL LS(A,N1,N2,NLINE,NCOL,B,ANS,ERR,ILL1,ILL2,AA,AB,AB2,AR)
*
      e=ans(1)
*
      do 600 i=1,jjj
      k=i+j40-1
      a(i,1)=str2(k)
      a(i,2)=1.
      b(i)=-fl(k)/area
*      write(6,*)a(i,1),a(i,2),b(i)
  600 continue
*
      CALL LS(A,N1,N2,NLINE,NCOL,B,ANS,ERR,ILL1,ILL2,AA,AB,AB2,AR)
*
      ed=ans(1)
*
      v=-e/ed
*
      do 633 i=1,jjj
      k=i+j40-1
      a(i,1)=-disp(k)/flen
      a(i,2)=1.
      b(i)=-fl(k)/area
*      write(6,*)a(i,1),a(i,2),b(i)
  633 continue
*
      CALL LS(A,N1,N2,NLINE,NCOL,B,ANS,ERR,ILL1,ILL2,AA,AB,AB2,AR)
*
      edisp=ans(1)
*      write(6,*)edisp
      delta=1./(e*100.)
      ddelta=1./edisp
      coeff=ddelta-delta
*
      do 8222 i=1,j60-1
      strrr(i)=-str1(i)
 8222 continue
      strrr1=-str1(i)
      strrr2=-((disp(j60)+dsta)/flen-fl(j60)/area*coeff)*100.
      dstrrr=strrr2-strrr1
      do 8223 i=j60,kend
      strrr(i)=-((disp(i)+dsta)/flen-fl(i)/area*coeff)*100.-dstrrr
 8223 continue
*
*      open(8,file='user')
      write(7,1000)smax,e*100/1000.,v,strrr(kmax),etc
 1000 format(1h ,'   qu(MPa)    E(GPa)         v    eac(%)    elc(%)',/,
     >       1h ,2f10.1,3f10.2)
      write(7,*)'pltファイルが作成されました'
*      write(8,*)'確認したら何かキーを押して下さい'
*      write(8,*)'続けて解析できます。'
*      iresult=incharqq()
*      close(8)
 2000 format(1h ,'ksta,dsta,coeff',/,1h ,i5,f10.3,g12.4)
      write(6,*)' rawファイルから変換したプロット用ファイル'
      write(6,*)
     > ' 時間（秒）、軸ひずみ（%）、周ひずみ（%）、軸応力（MPa）'
      do 1777 j=ksta-4,kend
*      write(*,*)j
      write(6,*)time(j),strrr(j),str2(j),-fl(j)/area
 1777 continue
      close(6)
      close(5)
      return
      end
*******************************************************
      SUBROUTINE LS(A,N1,N2,NLINE,NCOL,B,ANS,ERR,ILL1,ILL2,AA,AB,AB2,AR)
C   TO SOLVE 'A*ANS=B' USING LEAST-SQUARES METHOD
C   STANDARD DEVIATIONS ARE OUTPUTED TO 'ERR'
C
C  ILL1=10  NLINE IS LESS THAN NCOL
C       -1  NLINE EQUALS NCOL
C
      DIMENSION A(N1,N2),B(N1),ANS(N2),ERR(N2),
     +          AA(N2,N2),AB(N2),AB2(N2),AR(N2,N2)
*      WRITE(6,7676)B
 7676 FORMAT(1H ,'B IN LEAST',9G12.4)
      ILL1=0
      ILL2=0
      IF(NLINE.GE.NCOL)GOTO50
      ILL1=10
      WRITE(6,1000)ILL1
 1000 FORMAT(1H ,'ILL1=',I3,' NLINE IS LESS THAN NCOL',/,
     +1H ,'SOLUTION WAS NOT OBTAINED')
      RETURN
   50 DO 10 K=1,NCOL
      ERR(K)=0.
      AB(K)=0.
      DO 10 L=1,NCOL
      AA(K,L)=0.
   10 CONTINUE
      DO 20 K=1,NCOL
      DO 20 J=1,NLINE
      AB(K)=AB(K)+A(J,K)*B(J)
      DO 20 L=1,NCOL
      AA(K,L)=AA(K,L)+A(J,K)*A(J,L)
   20 CONTINUE
C     DO 5998 M=1,6
C     WRITE(6,5999)(AA(M,N),N=1,6)
C5999 FORMAT(1H ,'AA',6G12.4)
C5998 CONTINUE
C     DO 9998 M=1,10
C     WRITE(6,9999)(A(M,N),N=1,6)
C9999 FORMAT(1H ,'A',6G12.4)
C9998 CONTINUE
      DO 25 K=1,NCOL
      DO 25 J=1,NCOL
      AR(K,J)=AA(K,J)
   25 CONTINUE
      CALL SWEEP(AA,AB,AB2,N2,NCOL)
      ILL2=0
      IF(ILL2.EQ.0)GOTO30
      WRITE(6,1001)ILL2
 1001 FORMAT(1H ,'ILL2=',I3)
      RETURN
   30 DO 40 K=1,NCOL
      ANS(K)=AB2(K)
   40 CONTINUE
      IF(NLINE.GT.NCOL)GOTO60
 2005 ILL1=-1
      WRITE(6,1002)ILL1
 1002 FORMAT(1H ,'ILL1=',I3)
      WRITE(6,8111)
 8111 FORMAT(1H ,'ERROR IS NOT CALCULABLE')
      RETURN
C CALCULATE STANDARD DEVIATIONS
   60 DO 100 I=1,NLINE
      DO 100 K=1,NCOL
      B(I)=B(I)-A(I,K)*ANS(K)
  100 CONTINUE
*       WRITE(6,7474)(B(K),K=1,NLINE)
* 7474 FORMAT(1H ,'RESIDUAL',10G10.4)
 2003 IF(NLINE.EQ.NCOL)GOTO2005
      SVV=0.
      DO 105 I=1,NLINE
      SVV=SVV+B(I)**2
  105 CONTINUE
*       WRITE(6,2001)SVV
* 2001 FORMAT(1H ,'SVV=',G12.4)
      RVV=SQRT(SVV/FLOAT(NLINE-NCOL))
C
      DO 110 I=1,NCOL
      DO 120 J=1,NCOL
      IF(I.EQ.J)GOTO130
      AB(J)=0.
      GOTO120
  130 AB(J)=1.
  120 CONTINUE
      DO 140 K=1,NCOL
      DO 140 L=1,NCOL
  140 AA(K,L)=AR(K,L)
C
**      CALL DEB(3)
C
      CALL SWEEP(AA,AB,AB2,N2,NCOL)
C
*      CALL DEB(4)
C
      ILL2=0
      IF(ILL2.EQ.0)GOTO150
      WRITE(6,1001)ILL2
      RETURN
  150 ERR(I)=RVV*SQRT(ABS(AB2(I)))
  110 CONTINUE
      RETURN
      END
C
C      *****************************************************
C
      SUBROUTINE SWEEP(A,B,X,N2,N)
      DIMENSION A(N2,N2),B(N2),X(N2)
      NM1=N-1
      DO 1 K=1,NM1
      KP1=K+1
      DO 2 J=KP1,N
    2 A(K,J)=A(K,J)/A(K,K)
      B(K)=B(K)/A(K,K)
      DO 3 I=KP1,N
      DO 4 J=KP1,N
    4 A(I,J)=A(I,J)-A(I,K)*A(K,J)
    3 B(I)=B(I)-A(I,K)*B(K)
    1 CONTINUE
      X(N)=B(N)/A(N,N)
      DO 5 K=1,NM1
      L=N-K
      LP1=L+1
      X(L)=B(L)
      DO 6 J=LP1,N
    6 X(L)=X(L)-A(L,J)*X(J)
    5 CONTINUE
      RETURN
      END
*
* PRINCIPAL STRESS
*
      subroutine sig(status)
      use msflib
      logical status
      open(7,file='user')
      PAI=3.141593
      WRITE(7,*)'Input σx, σy, τxy'
      READ(7,*)SX,SY,SXY
        IF((SX-SY).EQ.0.)THEN
        IF(SXY.GT.0.)PH=PAI/4.
        IF(SXY.LT.0.)PH=-PAI/4.
        if(sxy.eq.0.)ph=0.
        ELSE
        PH=0.5*ATAN((2.*SXY)/(SX-SY))
          IF((SX-SY).LT.0.)PH=PH+PAI/2.
        ENDIF
      CP=COS(PH)
      SP=SIN(PH)
      S1=CP**2*SX+SP**2*SY+2.*SP*CP*SXY
      S2=SP**2*SX+CP**2*SY-2.*SP*CP*SXY
      PH=PH*180./PAI
      WRITE(7,100)S1,S2,PH
  100 FORMAT(1H ,'σ1,σ2,α',3g12.4)
*      write(7,*)'Hit any key'
*      iresult=incharqq()
*      close(7)
      return
      END
*
* coodinate conversion
*
      subroutine sig2(status)
      use msflib
      dimension s1(2,2),s2(2,2),c(2,2)
      logical status
      open(7,file='user')
      PAI=3.141593
      WRITE(7,*)'Input σx, σy, τxy, α'
      READ(7,*)SX,SY,SXY,alph
      alpha=alph*pai/180.
      sa=sin(alpha)
      ca=cos(alpha)
      s1(1,1)=sx
      s1(1,2)=sxy
      s1(2,1)=sxy
      s1(2,2)=sy
      c(1,1)= ca
      c(1,2)= sa
      c(2,1)=-sa
      c(2,2)= ca
      do 100 i=1,2
      do 100 j=1,2
      s2(i,j)=0.
  100 continue
      do 200 i=1,2
      do 200 j=1,2
      do 200 k=1,2
      do 200 l=1,2
      s2(i,l)=s2(i,l)+c(i,j)*s1(j,k)*c(l,k)
  200 continue
      WRITE(7,300)s2(1,1),s2(2,2),s2(1,2)
  300 FORMAT(1H ,'σx'',σy'',τx''y''',3g12.4)
*      write(7,*)'Hit any key'
*      iresult=incharqq()
*      close(7)
      return
      END
*
* Principal stress in three dimension
*
      subroutine prst(status)
      use msflib
      dimension sig(6),pst(3),dir(3,2),dcc(3,3)
      logical status
      open(7,file='user')
      pai=3.141593
      WRITE(7,4504)
      write(7,*)'Input σx,σy,σz,τxy,τyz,τzx'
      read(7,*)(sig(i),i=1,6)
      CALL FPRST(SIG,PST,DIR,DCC)
 4504 FORMAT(' Principal stresses and direction(degree),',/,
     >       ' phai is measured from x to y',/,
     >       ' theta is measured from z to x-y plane',/,
     >       ' To draw stereo projection (lower hemisphere),',/,
     >       ' draw a circle whose radius is 2.5cm',/,
     >       ' and plot points given below as (x, y)')
      WRITE(7,4507)
      DO 3801 KKk=1,3
      kk=4-kkk
      PHAI=DIR(KK,1)*PAI/180.
      THETA=DIR(KK,2)*PAI/180.
      R=2.5*TAN(THETA/2.)
      X=-R*cos(PHAI)
      Y=-R*sin(PHAI)
      WRITE(7,3800)KKk,PST(KK),DIR(KK,1),DIR(KK,2),X,Y
 3801 CONTINUE
 4507 FORMAT('         STRESS      PHAI     THETA         X         Y')
 3800 FORMAT(' SIG',I1,f10.4,4F10.1)
*
*      write(7,*)'Hit any key'
*      iresult=incharqq()
*      close(7)
      return
      END
C*****************************************************************
      SUBROUTINE FPRST(BB,PST,DIR,DCC)
C  BB IS SIX-INDEPENDENT COMPONENTS OF SECOND RANKED
C  TENSOR (XX,YY,ZZ,XY,YZ,ZX).
C  PST IS ITS PRINCIPAL VALUE, DIR IS ITS ANGLE, DCC IS
C  ITS DIRECTION COSINE.
C   PST(K),DIR(K,J),DCC(K,L)
C   K=1:MIN 2:MID 3:MAX
C   J=1:PHAI (FROM X TO Y) 2:THETA (AGAINST Z)
C   L=1:AGAINST X 2:AGAINST Y 3:AGAINST Z
      DIMENSION A(3,3),V(3,3),E(3),BB(6),PST(3),DIR(3,2),DE(3,2)
     + ,DCC(3,3),DC(3,3)
      A(1,1)=BB(1)
      A(1,2)=BB(4)
      A(1,3)=BB(6)
      A(2,1)=BB(4)
      A(2,2)=BB(2)
      A(2,3)=BB(5)
      A(3,1)=BB(6)
      A(3,2)=BB(5)
      A(3,3)=BB(3)
      CALL FEV(A,E,V,3,100)
      DO 200 K=1,3
      FL=V(1,K)
      FM=V(2,K)
      FN=V(3,K)
      FKK=FL*FL+FM*FM+FN*FN
      FK2=SQRT(FKK)
      FL=FL/FK2
      FM=FM/FK2
      FN=FN/FK2
      DC(K,1)=FL
      DC(K,2)=FM
      DC(K,3)=FN
      CALL FANG(FL,FM,FN,PHAI,THETA)
      DE(K,1)=PHAI
      DE(K,2)=THETA
  200 CONTINUE
      DO 345 K=1,3
      IF(DC(K,3).GE.0)GOTO345
      DO 346 J=1,3
      DC(K,J)=-DC(K,J)
  346 CONTINUE
  345 CONTINUE
      AMIN=1.E30
      DO 3 K=1,3
      IF(E(K).GT.AMIN)GOTO3
      AMIN=E(K)
      KMIN=K
    3 CONTINUE
      PST(1)=E(KMIN)
      DIR(1,1)=DE(KMIN,1)
      DIR(1,2)=DE(KMIN,2)
      DCC(1,1)=DC(KMIN,1)
      DCC(1,2)=DC(KMIN,2)
      DCC(1,3)=DC(KMIN,3)
      KN=0
      DO 4 K=1,3
      IF(K.EQ.KMIN)GOTO4
      KN=KN+1
      E(KN)=E(K)
      DE(KN,1)=DE(K,1)
      DE(KN,2)=DE(K,2)
      DC(KN,1)=DC(K,1)
      DC(KN,2)=DC(K,2)
      DC(KN,3)=DC(K,3)
    4 CONTINUE
      IF(E(1).LT.E(2))GOTO5
      PST(2)=E(2)
      DIR(2,1)=DE(2,1)
      DIR(2,2)=DE(2,2)
      DCC(2,1)=DC(2,1)
      DCC(2,2)=DC(2,2)
      DCC(2,3)=DC(2,3)
      PST(3)=E(1)
      DIR(3,1)=DE(1,1)
      DIR(3,2)=DE(1,2)
      DCC(3,1)=DC(1,1)
      DCC(3,2)=DC(1,2)
      DCC(3,3)=DC(1,3)
      GOTO6
    5 PST(2)=E(1)
      DIR(2,1)=DE(1,1)
      DIR(2,2)=DE(1,2)
      DCC(2,1)=DC(1,1)
      DCC(2,2)=DC(1,2)
      DCC(2,3)=DC(1,3)
      PST(3)=E(2)
      DIR(3,1)=DE(2,1)
      DIR(3,2)=DE(2,2)
      DCC(3,1)=DC(2,1)
      DCC(3,2)=DC(2,2)
      DCC(3,3)=DC(2,3)
    6 CONTINUE
      RETURN
      END
C********************************************
      SUBROUTINE FANG(FL,FM,FN,PHAI,THETA)
      PAI2=1.57080
      RTD=57.2958
      THR=0.000001
      IF(FN.GT.0.)GOTO10
      FN=-FN
      FL=-FL
      FM=-FM
   10 AFN=ABS(FN)
      AFL=ABS(FL)
      AFM=ABS(FM)
      IF(AFL.GT.THR)GOTO20
      PHAI=PAI2
      IF(AFN.GT.THR)GOTO30
      THETA=PAI2
      GOTO1000
   30 THETA=ATAN(FM/FN)
      GOTO1000
C
   20 IF(AFM.GT.THR)GOTO120
      PHAI=0.
      IF(AFN.GT.THR)GOTO130
      THETA=PAI2
      GOTO1000
  130 THETA=ATAN(FL/FN)
      GOTO1000
C
  120 IF(AFN.GT.THR)GOTO220
      THETA=PAI2
      IF(AFL.GT.THR)GOTO230
      PHAI=PAI2
      GOTO1000
  230 PHAI=ATAN(FM/FL)
      GOTO1000
C
C
  220 PHAI=ATAN(FM/FL)
      THETA=ATAN(FL/(FN*COS(PHAI)))
 1000 IF(THETA.GE.0)GOTO6111
      PHAI=PHAI+PAI2*2.
      THETA=-THETA
 6111 PHAI=PHAI*RTD
      THETA=THETA*RTD
      RETURN
      END
C **********************************************
      SUBROUTINE FEV(A,E,V,N,KMAX)
      DIMENSION A(3,3),V(3,3),E(3)
      INTEGER P,Q
      NM1=N-1
      DO 1 I=1,N
      DO 2 J=1,N
    2 V(I,J)=0.
    1 V(I,I)=1.
      DO 100 KAISU=1,KMAX
      AMAX=0.
      DO 5 I=1,NM1
      IP1=I+1
      DO 5 J=IP1,N
      IF(ABS(A(I,J)).LT.AMAX)GOTO5
      P=I
      Q=J
      AMAX=ABS(A(P,Q))
    5 CONTINUE
      IF(ABS(A(P,Q)).LT.1.0E-05)GOTO73
      T=0.785398
      BUNBO=A(P,P)-A(Q,Q)
      IF(ABS(BUNBO).LT.1.0E-10)GOTO10
      R=2.0*A(P,Q)/BUNBO
      T=0.5*ATAN(R)
   10 S=SIN(T)
      C=COS(T)
      DO 11 J=1,N
      APJ=A(P,J)
      AQJ=A(Q,J)
      A(P,J)=APJ*C+AQJ*S
   11 A(Q,J)=-APJ*S+AQJ*C
      DO 12 I=1,N
      AIP=A(I,P)
      AIQ=A(I,Q)
      A(I,P)=AIP*C+AIQ*S
      A(I,Q)=-AIP*S+AIQ*C
      VIP=V(I,P)
      VIQ=V(I,Q)
      V(I,P)= VIP*C+VIQ*S
   12 V(I,Q)=-VIP*S+VIQ*C
  100 CONTINUE
   73 DO 13 I=1,N
   13 E(I)=A(I,I)
      RETURN
      END

