      program temp
      
C==注意　この計算プログラムの結果は計算精度的に問題がある場合があります=
C==論文等に利用する場合は適切に書き換えてください=======================

C==プログラム製作者である上田はプログラミングに関しては全くの素人です===
C==プログラム製作者である上田はプログラミングに関しては全くの素人です===
C==プログラム製作者である上田はプログラミングに関しては全くの素人です===
C==必要があれば適切に書き換えてください=================================
C=======================================================================
C=======================================================================
C=========使い方========================================================
C42-47行目の条件を入力する
C55,56行目の圧力と流速を複素数で入力する．
      IMPLICIT NONE
      integer b,bmax
      DOUBLE PRECISION cp,rhom,Mu,Nyu,Kappa,Alpha,Pr,Ss,Tm,Gam,Pm
      DOUBLE PRECISION R,PI,det,f
      DOUBLE PRECISION Tc,TH,dTm
      DOUBLE PRECISION Lstack,DelX,w
      DOUBLE PRECISION wtu,wta
      COMPLEX*16 Pin,Uin,Pout,Uout,P,U
      DOUBLE PRECISION Iin,Iout,In
      COMPLEX*16 Ya,Yu,Xa,Xu,I  
      COMPLEX*16 ME(2,2),Ms(2,2), Mf(2,2),Mff(2,2),Msta(2,2)
      COMPLEX*16 DJ1J0      
      
C=================================================
      OPEN(54,FILE='G.TXT')
      bmax=100
C===========PiとIの定義=====================================
      PI=4D0*DATAN(1D0)      !円周率の定義
      I =CMPLX(0D0,1D0)      !虚数項の定義
C===========単位行列の定義==================================
      ME(1,1)=CMPLX(1E0,0E0)
      ME(1,2)=CMPLX(0E0,0E0)
      ME(2,1)=CMPLX(0E0,0E0)
      ME(2,2)=CMPLX(1E0,0E0)  

C=====================実験条件====================
C===========================================================
      R=0.34D-3   !蓄熱器の流路半径R粘性                                 
      Pm=1*101E3         !平均気圧
      Tc=293          !蓄熱器低温端の温度
      TH=293*2.3        !蓄熱器高温端の温度
      Lstack=20E-3      !蓄熱器の長さ
      f=100D0             !周波数
c====================入り口の計算====================================== 

      w=2D0*f*Pi
      
      Tm=TC
      CALL bussei(cp,rhom,Mu,Nyu,Kappa,Alpha,Pr,Ss,Tm,Gam,Pm)
      
      Pin=cmplx(1.08371789371,-0.19130178936)     
      Uin=cmplx(1.777853001082D-03,8.733481219974D-05)                                                 
      
      Iin=0.5*DBLE(CONJG(Pin)*Uin)

      Mff(1,1)=ME(1,1)
      Mff(2,1)=ME(2,1)
      Mff(1,2)=ME(1,2)
      Mff(2,2)=ME(2,2)      
c==========================================================
      DelX=Lstack/bmax        !管の刻み幅
      dTm=(TH-TC)/Lstack      !温度勾配
      
      DO 3 b=1,bmax           !蓄熱器内の計算スタート
      Tm=TC+dTm*DelX*b
      CALL bussei(cp,rhom,Mu,Nyu,Kappa,Alpha,Pr,Ss,Tm,Gam,Pm) !物性値の計算


      wtu =R**2*w/(2D0*Nyu )              !このあたりの記号の使い方はは富永著の熱音響工学の基礎を参照
      wta =R**2*w/(2D0*Alpha )
      
      Ya =CMPLX(1D0,1D0)*CMPLX(DSQRT(wta),0D0)
      Yu =CMPLX(1D0,1D0)*CMPLX(DSQRT(wtu),0D0)            
      
      Xa =2D0/I/Ya*DJ1J0(I*Ya)
      Xu =2D0/I/Yu*DJ1J0(I*Yu)

      Ms(1,1)=0D0
      Ms(1,2)=(-I*w*rhom)/(1D0-Xu)
      Ms(2,1)=-I*w*(1D0+(Gam-1D0)*Xa)/(Gam*Pm)
      Ms(2,2)=(Xa-Xu)*dTm/(1D0-Xu)/(1D0-Pr)/Tm

      Mf(1,1)=ME(1,1)+delX*Ms(1,1)
      Mf(1,2)=ME(1,2)+delX*Ms(1,2)
      Mf(2,1)=ME(2,1)+delX*Ms(2,1)
      Mf(2,2)=ME(2,2)+delX*Ms(2,2)
      

      CALL MATBY(Mf,Mff,Msta)


      Mff(1,1)=Msta(1,1)
      Mff(2,1)=Msta(2,1)
      Mff(1,2)=Msta(1,2)
      Mff(2,2)=Msta(2,2)
      P=Pin*Mff(1,1)+Uin*Mff(1,2)
      U=Pin*Mff(2,1)+Uin*Mff(2,2)
      In=0.5*DBLE(P*CONJG(U))
      write(54,*)CDABS(P),CDABS(U),In/Iin
3     continue
      
      
           
      Pout=Mff(1,1)*Pin+Mff(1,2)*Uin
      Uout=Mff(2,1)*Pin+Mff(2,2)*Uin
      Iout=0.5*DBLE(CONJG(Pout)*Uout)
      
      write(*,*)'増幅率',Iout/Iin,'周波数',f
      write(*,*)'出口の音響インピーダンス',Pout/Uout
      write(*,*)'出口での特性音響インピーダンス',rhom*Ss

      stop
      End
      
      
C=========物性値の計算======================================
      SUBROUTINE bussei(cp,rhom,Mu,Nyu,Kappa,Alpha,Pr,Ss,Tm,Gam,Pm)
      
      DOUBLE PRECISION cp,rhom,Mu,Nyu,Kappa,Alpha,Pr,Ss,Tm,Gam,Pm
      
           

      
      
      rhom=(4.1461D0
     &       -0.019705*Tm
     &       +4.7881e-05*Tm*Tm
     &       -6.2842e-08*Tm*Tm*Tm
     &       +4.2339e-11*Tm*Tm*Tm*Tm
     &       -1.1472e-14*Tm*Tm*Tm*Tm*Tm)*Pm/101D3

      Mu=-7.0031D-06
     &   +0.14018D-06*Tm
     &   -0.00026899D-6*Tm*Tm
     &   +3.5727D-13*Tm*Tm*Tm
     &   -2.4676D-16*Tm*Tm*Tm*Tm
     &   +6.8396D-20*Tm*Tm*Tm*Tm*Tm  
     
     
      Kappa=4.9135D-3
     &       +0.053336D-3*Tm
     &       +0.00013329D-3*Tm*Tm
     &       -3.4467D-10*Tm*Tm*Tm
     &       +3.4584D-13*Tm*Tm*Tm*Tm
     &       -1.2551D-16*Tm*Tm*Tm*Tm*Tm
     
      Gam=1.3803D0
     &      +0.00020301*Tm
     &      -5.1868D-7*Tm*Tm
     &      +2.4927D-10*Tm*Tm*Tm
     &      +1.0045D-13*Tm*Tm*Tm*Tm
     &      -7.8378D-17*Tm*Tm*Tm*Tm*Tm
C      WRITE(*,*)Tm,Gam
      Cp=1057.9D0
     &      -0.3684*Tm
     &      +0.00063128*Tm*Tm
     &      +3.4653D-7*Tm*Tm*Tm
     &      -8.7851D-10*Tm*Tm*Tm*Tm 
     &      +3.6441D-13*Tm*Tm*Tm*Tm*Tm 
C      WRITE(*,*)Tm,rhom,Mu,Kappa     
      Nyu =Mu /rhom
      Alpha = Kappa / rhom / Cp
      Pr =Mu / Kappa *Cp
      Ss =DSQRT(Gam * Pm /rhom ) 
c      WRITE(*,*)Tm,Pr
      RETURN
      END
      

     
     
     
C*******************************************************************************
C=========行列の掛け算======================================
      SUBROUTINE MATBY(Maa,Mbb,Mcc)
      COMPLEX*16 Maa(2,2),Mbb(2,2),Mcc(2,2)
      
      Mcc(1,1)=Maa(1,1)*Mbb(1,1)+Maa(1,2)*Mbb(2,1)
      Mcc(1,2)=Maa(1,1)*Mbb(1,2)+Maa(1,2)*Mbb(2,2)
      Mcc(2,1)=Maa(2,1)*Mbb(1,1)+Maa(2,2)*Mbb(2,1)
      Mcc(2,2)=Maa(2,1)*Mbb(1,2)+Maa(2,2)*Mbb(2,2)
      
      RETURN
      END
C==============================================================


C*******************************************************************************
C     FUNCTION OF J1(Z)/J0(Z) DEVELOPED BY SEIICHI WASHIO
C     J1(Z), J0(Z): 1ST AND 0-TH ORDER BESSEL FUNCTIONS OF THE FIRST KIND
C     DEVELOPED BY SEIICHI WASHIO
C*******************************************************************************
C
      COMPLEX*16 FUNCTION DJ1J0(Z)
C      
      IMPLICIT COMPLEX*16 (E,F,Z)
      IMPLICIT DOUBLE PRECISION (A,D,X,Y) 
      COMPLEX*16 Z

      IF(CDABS(Z).GT.26.0) GO TO 9
      M=1.581*CDABS(Z)+1.0
      Z1=Z/2.0D+00
      Z2=Z1*Z1
      J1=0
      J0=0
      K=0
      E1=1.0D+00
      E0=1.0D+00
      F1=1.0D+00
      F0=1.0D+00
      AK1=1.0D+00
    1 K=K+1
      AK=AK1
      AK1=AK1+1.0D+00
      F0=-F1*Z2/AK
      F1=F0/AK1
      IF(K.GT.M)GO TO 2
      E1=E1+F1
      E0=E0+F0
      GO TO 1
    2 IF(J1.EQ.1)GO TO 5
      S1=CDABS(F1)
      IF(ABS(S1/SNGL(DREAL(E1))).GE.1.0D-16)GO TO 3
      IF(ABS(S1/SNGL(DIMAG(E1))).LT.1.0D-16)GO TO 4
    3 E1=E1+F1
      IF(J0.EQ.1)GO TO 1
      GO TO 5
    4 J1=1
      IF(J0.EQ.1)GO TO 8
    5 S0=CDABS(F0)
      IF(ABS(S0/SNGL(DREAL(E0))).GE.1.0D-16)GO TO 6
      IF(ABS(S0/SNGL(DIMAG(E0))).LT.1.0D-16)GO TO 7
    6 E0=E0+F0
      GO TO 1
    7 J0=1
      IF(J1.EQ.1)GO TO 8
      GO TO 1
    8 DJ1J0=Z1*E1/E0
      RETURN
C     APPROXIMATION BY HANKEL'S ASYMPTOTIC EXPANSION
    9 X=DREAL(Z)
      Y=DIMAG(Z)
      I=1 
      IF(X.GE.0.0)GO TO 10 
      X=-X
      Y=-Y
      I=2
   10 Z1=8.0D+00*DCMPLX(X,Y)
      ZE0=1.0D+00
      E0=1.0D+00
      ZE1=1.0D+00
      E1=1.0D+00
      ZF0=-1.0D+00/Z1
      F0=-1.0D+00/Z1
      ZF1=3.0D+00/Z1
      F1=3.0D+00/Z1
      Z2=Z1**2
      IF(Y.LT.0.0)GO TO 11
      D1=DEXP(-2.0D+00*Y)
      D0=1.0D+0-D1
      D1=1.0D+0+D1
      GO TO 12
   11 D1=DEXP(2.0D+00*Y)
      D0=D1-1.0D+00
      D1=1.0D+00+D1
   12 AC=DCOS(X)
      AS=DSIN(X)
      DO 13 L=1,8
      A1=DFLOAT((4*L-3)**2)
      A2=DFLOAT((4*L-1)**2)
      A3=DFLOAT(2*L*(2*L-1))
      ZE0=-ZE0*A1*A2/Z2/A3
      ZE1=-ZE1*(4.0D+00-A1)*(4.0D+00-A2)/Z2/A3
      E0=E0+ZE0
   13 E1=E1+ZE1
      DO 14 L=1,7
      A1=DFLOAT((4*L-1)**2)
      A2=DFLOAT((4*L+1)**2)
      A3=DFLOAT(2*L*(2*L+1))
      ZF0=-ZF0*A1*A2/Z2/A3
      ZF1=-ZF1*(4.0D+00-A1)*(4.0D+00-A2)/Z2/A3
      F0=F0+ZF0
   14 F1=F1+ZF1
      ZA=DCMPLX(D1*AC,-D0*AS)
      ZB=DCMPLX(D1*AS,D0*AC)
      DJ1J0=((F1-E1)*ZA+(E1+F1)*ZB)/((E0+F0)*ZA+(E0-F0)*ZB)
      IF(I.EQ.1)RETURN
      DJ1J0=-DJ1J0
      RETURN
      END
     
