
program streamf

  use dcl
  real :: k,U,mu,Ci,Cr,c,z
  integer,parameter :: nx=50, ny=50
  real,parameter :: xmin=0., xmax=2.5*DCL_PI, ymin=0., ymax=1.
  real,dimension(0:nx,0:ny) :: psi
  
  U=1.0
!  mu=1.61
  mu=2.5
  Cr=U*0.5
!  Ci=0.31*U/mu
  Ci=0	
  c2=Cr**2+Ci**2
  
!-- 球面調和関数 ----
open(17, file='temp.dat', status='unknown')
    do j=0,ny
       do i=0,nx
        k = ( xmin + (xmax-xmin)*i/nx )
        z = ( ymin + (ymax-ymin)*j/ny )
!        if (i==0) write (*,*) z
        psi(i,j) = mu*cos(k)*(sinh(mu*z)-U*Cr*cosh(mu*z)/(mu*c2))&
             -mu*sin(k)*U*Ci*cosh(mu*z)/(mu*c2)
        write (17,'(3e12.4)') k, z, psi(i,j) 
      end do
    end do
close(17)

!-- グラフ ----
    call swcset('FNAME','hoge')
    call DclOpenGraphics()
    call DclNewFrame

    call DclSetWindow( xmin, xmax, ymin, ymax )
    call DclSetViewPort( 0.2, 1.0, 0.2, 0.5 )
    call DclSetTransFunction

    call DclDrawScaledAxis

    call DclDrawContour( psi )

    call DclCloseGraphics

end program
