C使い方
C34-41行目の装置の条件に値を入れる
C46-47行目の境界条件を入れる
Cコンパイル実行
Cacoustic_field.txt と work_flow.txtが生成 前者には位置，圧力振幅，流速振幅，圧力の位相，後者には位置，仕事流が書き込まれる
C5    プログラムを書くことに関してはど素人のため，スマートな命令分ではない





C=====================================================================================================      
      program tube
      IMPLICIT NONE
      integer N,Nmax
      COMPLEX*16 P0,U0,P,U
      DOUBLE PRECISION Pamp,Uamp,Phi,W,Pphi,Uphi      
      COMPLEX*16 I,Ya,Yu,Xa,Xu,k,Mtube(2,2),Mtube1(2,2)
      DOUBLE PRECISION Pi,r,L,A,TH,k0,omega,f,cp,rhom
      DOUBLE PRECISION Mu,Nyu,Kappa,Alpha,Pr,Ss,Tm,Gam,Pm
      DOUBLE PRECISION wta,wtu,X,dx
      COMPLEX*16 DJ1J0
C============================================================
      OPEN(10,FILE='acoustic_field.TXT')
      OPEN(11,FILE='work_flow.TXT')
C============================================================
      Pi=4D0*DATAN(1D0)       !円周率の定義
      I =CMPLX(0D0,1D0)       !複素数iの定義
      Nmax=100              !管の分割数
C============================================================
C=========管の形状など=======================================
C============================================================
      Pm=101D3        !mean pressure
      r=10.5D-3         !radius of tube
      L=1.00D0          !length of tube
      A=r*r*Pi          !cross-sectional area of tube
      TH=300D0          !mean temperature
      f=148.5D0         !frequency
      omega=2D0*Pi*f    !angular frequency
      dx=-L/real(Nmax)  !ここの-をなくせば，正の方向に進む．***重要事項***

C============================================================
C=========境界条件=======================================
C============================================================   
      P0=CMPLX(0.78D3,0D0)     !pressure at tube-end 閉端の場合は有限値，開端の場合はゼロ
      U0=CMPLX(0D0,0D0)       !velocity at tube-end 開端の場合はゼロ  , 開端の場合は有限値
      Tm=TH
      CALL bussei(cp,rhom,Mu,Nyu,Kappa,Alpha,Pr,Ss,Tm,Gam,Pm)    !物性値の読み込み，初期設定は空気

                 
      wta =(r*r)*omega/(2D0*Alpha )         !おめが　たう　アルファ
      wtu =(r*r)*omega/(2D0*Nyu )           !おめが　たう　ニュー
  
      Ya =CMPLX(1D0,1D0)*CMPLX(DSQRT(wta),0D0)
      Yu =CMPLX(1D0,1D0)*CMPLX(DSQRT(wtu),0D0)            

      Xa =2D0/I/Ya*DJ1J0(I*Ya)            !かい　アルファー　DJ1J0はベッセル関数の比
      Xu =2D0/I/Yu*DJ1J0(I*Yu)            !かい　ニュー　DJ1J0はベッセル関数の比

      K0=omega/Ss                                     !自由空間中の波数
      K=k0*CSQRT( ( 1D0+(Gam-1D0)*Xa )  /(1D0-Xu) )   !半径rの管中の波数
      
      Do 100 N=1,Nmax          !長さLの管をNmax分割して音場を計算
      X=(real(N)*dx)
      
      Mtube(1,1)=CDCOS(K*X)                           !閉端からXの位置までの管の伝達行列1 1成分
      Mtube(1,2)=-(I*omega*rhom)/(K*(1-Xu))*CDSIN(K*X) !閉端からXの位置までの管の伝達行列1 2成分
      Mtube(2,1)=(K*(1-Xu))/(i*omega*rhom)*CDSIN(K*X)  !閉端からXの位置までの管の伝達行列2 1成分
      Mtube(2,2)=CDCOS(K*X)                           !閉端からXの位置までの管の伝達行列2 2成分
      
      P=P0*Mtube(1,1)+U0*Mtube(1,2)       !位置xでの圧力(複素数）
      U=P0*Mtube(2,1)+U0*Mtube(2,2)       !位置xでの流速(複素数)
      Pamp=CDABS(P)                       !位置xでの圧力振幅(実数)
      Uamp=CDABS(U)                       !位置xでの流速振幅(実数)
      
      Pphi=DATAN2(AIMAG(P),DBLE(P))       !位置xでの圧力の位相(実数，ラジアン)，基準はP0
      Uphi=DATAN2(AIMAG(U),DBLE(U))       !位置xでの流速の位相(実数，ラジアン)，基準はP0
      Phi=(Uphi-Pphi)                     !位置xでの圧力から見た流速の位相（実数，ラジアン）
                           
      W=0.5*DBLE(P*CONJG(U))              !位置xでの仕事流，断面積Aを掛ければpower flowになる.
      write(10,*)x,Pamp,Uamp,Pphi         !ファイル10（acoustic_field.txtに（位置，圧力振幅，流速振幅，圧力の位相）を書きこみ
      write(11,*)x,W                      !ファイル10（work_flow.txtに(位置，仕事流(W/m^2))を書きこみ

100   CONTINUE
      close(10)
      close(11)
      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
      
           

      
C密度
      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
C粘性係数
      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  
     
C熱伝導度
      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
C比熱比
      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定圧比熱
      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動粘性係数     
      Nyu =Mu /rhom
C熱拡散係数
      Alpha = Kappa / rhom / Cp
Cプラントル数
      Pr =Mu / Kappa *Cp
C音速
      Ss =DSQRT(Gam * Pm /rhom ) 
      RETURN
      END

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     このFunctionはベッセル関数の比を計算するのに使います．少々，tanhを使用する場合（平行平板）に比べて計算負荷が増える
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