
program main
    use SeismicAnalysisClass
    use LoggerClass
    implicit none

    type(SeismicAnalysis_)  :: siml 

    type(FEMDomain_),target :: domains(1)
    real(real64),allocatable :: x_axis(:),y_axis(:),z_axis(:)
    logical :: connected(1,1) = .false. 

    real(real64) ::  &
        Length, Width, Thickness
    real(real64),allocatable ::YoungModulus(:),PoissonRatio(:),rho(:)
    real(real64),allocatable :: v(:,:)
    
    real(real64) :: dt,t,angle
    type(MPI_) :: mpid
    
    real(real64) :: Density(1:4),Vs(1:4),Vp(1:4),H(1:4,1:2),E,vv
    integer(int32) :: layer
    type(Random_) :: random
    type(Logger_),allocatable :: loggers(:)
    type(IO_) :: file
    call mpid%start()

    call file%open("noize.txt")
    do i_i=1,10000
        call file%write( dble(i_i)/1000.0d0,random%gauss(mu=0.0d0,sigma=1.0d0) ) 
    enddo
    call file%close()
    stop
    ! 常時微動アレイ観測の再現シミュレーション

    ! 地盤データ
    ! rho (t/m^3)    Vp (m/s)    Vs (m/s)     Thickness (m)
    ! 1.4   1000    170   8.0
    ! 1.5   1100    250   4.0
    ! 1.65  1400    350   3.0
    ! 2.00  1900    600  \infty 

    ! 速度構造の入力
    density(1) =1.40d0; Vp(1) = 1000.0d0 ; Vs(1) = 170.0d0 ; H(1,1:2) = [      -8.0d0,  0.0d0 ];
    density(2) =1.50d0; Vp(2) = 1100.0d0 ; Vs(2) = 250.0d0 ; H(2,1:2) = [H(1,1)-4.0d0,  H(1,1)];
    density(3) =1.65d0; Vp(3) = 1400.0d0 ; Vs(3) = 350.0d0 ; H(3,1:2) = [H(2,1)-3.0d0,  H(2,1)];
    density(4) =2.00d0; Vp(4) = 1900.0d0 ; Vs(4) = 600.0d0 ; H(4,1:2) = [H(3,1)-3.0d0,  H(3,1)];

    Length    = 300.00d0 ! 1 km
    Width     = 300.00d0 ! 1 km
    Thickness = 30.00d0!  300 m
    !call domains(1)%create("Cube3D",x_num=300,y_num=300,z_num=30)
    call domains(1)%create("Cube3D",x_num=100,y_num=100,z_num=30)
    call domains(1)%resize(x=Length,y=width,z=Thickness)
    call domains(1)%move(z = - domains(1)%z_max() )

    call domains(1)%vtk("mesh")

    print *, domains(1)%nn()
    
    ! ヤング率およびポワソン比の計算
    YoungModulus = zeros(domains(1)%ne())
    PoissonRatio = zeros(domains(1)%ne())
    Rho          = zeros(domains(1)%ne())

    do layer=1,4
        vv = (((Vp(layer)/Vs(layer))**2) - 2.0d0)/2.0d0/((((Vp(layer)/Vs(layer))**2) - 1.0d0))
        E = (Vs(layer)*Vs(layer)*density(layer) )*2.0d0*(1.0d0+vv)
        print *, "LAYER " + str(layer) + " E(GPa) :",(E/1000.0d0/1000.0d0), " v :",vv,"Height: ",H(layer,2)
        YoungModulus(domains(1)%getElementList(zmax=H(layer,2)) ) = E 
        PoissonRatio(domains(1)%getElementList(zmax=H(layer,2)) ) = vv
        Rho(domains(1)%getElementList(zmax=H(layer,2)) ) = density(layer)
    enddo
    call domains(1)%vtk("YoungModulus",scalar=YoungModulus)
    call domains(1)%vtk("PoissonRatio",scalar=PoissonRatio)
    call domains(1)%vtk("Density",scalar=Rho)
    
    ! mesh creation done!
    
    
    ! >> Run solver
    print *, "Initialize >> "
    call siml%init(domains)

    print *, "Material >> "
    call siml%setMaterial(DomainID=1,density=Rho,YoungModulus=YoungModulus,PoissonRatio=PoissonRatio)
    
    print *, "Boundary >> "
    call siml%setBoundary(DomainID=1,NodeList=domains(1)%getNodeList(&
        zmax = domains(1)%zmin() ),condition="Absorbing Boundary", boundaryValue=[100.0d0, 1.00d0])
    call siml%setBoundary(DomainID=1,NodeList=domains(1)%getNodeList(&
        xmax = domains(1)%xmin() ),condition="Absorbing Boundary", boundaryValue=[100.0d0, 1.00d0])
    call siml%setBoundary(DomainID=1,NodeList=domains(1)%getNodeList(&
        ymax = domains(1)%ymin() ),condition="Absorbing Boundary", boundaryValue=[100.0d0, 1.00d0])
    call siml%setBoundary(DomainID=1,NodeList=domains(1)%getNodeList(&
        xmin = domains(1)%xmax() ),condition="Absorbing Boundary", boundaryValue=[100.0d0, 1.00d0])
    call siml%setBoundary(DomainID=1,NodeList=domains(1)%getNodeList(&
        ymin = domains(1)%ymax() ),condition="Absorbing Boundary", boundaryValue=[100.0d0, 1.00d0])

    print *, "Solve >> "

    siml%modal%solver%debug= .true.
    !siml%modal%solver%er0  = dble(1.0e-15)
    !siml%modal%solver%relative_er  = dble(1.0e-10)
    ! 1 kHz sampling
    dt = 1.0d0/1000.0d0

    allocate(loggers(37) )
    v = zeros(domains(1)%nn(),domains(1)%nd() )
    do j_j=1,36
        call Loggers(j_j)%set(&
            femdomain=domains(1),&
            position=[ 150.0d0+cos(radian(10.0d0*j_j) ),150.0d0+sin(radian(10.0d0*j_j)),0.0d0 ],&
            dataset=siml%A,&
            name="Array_A_"+zfill(j_j,4) )
        call loggers(j_j)%vtk("Array_A_"+zfill(j_j,4) )
        
    enddo
    j_j = 37
    call Loggers(j_j)%set(&
        femdomain=domains(1),&
        position=[ 150.0d0,150.0d0,0.0d0 ],&
        dataset=siml%A,&
        name="Array_A_"+zfill(j_j,4) )
    call loggers(j_j)%vtk("Array_A_"+zfill(j_j,4) )

    do j_j=1,37
        call loggers(j_j)%start()
    enddo
    ! damping :: default
    !siml%alpha = 0.0d0
    !siml%beta = 0.0d0
    
    t = 0.0d0
    do i_i=1,10000
        
        print *, t

        ! m/s
        ! 加速度，ホワイトノイズで遠方(0,0,0),(300,0,0),(0,300,0),(300,3000,0)を上下に揺らす.
        ! σ=1,平均0,倍率50.0d0/1000.0d0 (m/s/s)
        call siml%setBoundary(DomainID=1,NodeList=domains(1)%getNodeList(&
            xmax=1.0d0,ymax=1.0d0,zmin=-0.0010d0)&
            ,condition="A",boundaryValue=[0.0d0,0.0d0,500.0d0/1000.0d0*random%gauss(mu=0.0d0,sigma=1.0d0) ])
        
        call siml%setBoundary(DomainID=1,NodeList=domains(1)%getNodeList(&
            xmin=299.0d0,ymax=1.0d0,zmin=-0.0010d0)&
            ,condition="A",boundaryValue=[0.0d0,0.0d0,500.0d0/1000.0d0*random%gauss(mu=0.0d0,sigma=1.0d0) ])
        
        call siml%setBoundary(DomainID=1,NodeList=domains(1)%getNodeList(&
            xmax=1.0d0,ymin=299.0d0,zmin=-0.0010d0)&
            ,condition="A",boundaryValue=[0.0d0,0.0d0,500.0d0/1000.0d0*random%gauss(mu=0.0d0,sigma=1.0d0) ])
        
        call siml%setBoundary(DomainID=1,NodeList=domains(1)%getNodeList(&
            xmin=299.0d0,ymin=299.0d0,zmin=-0.0010d0)&
            ,condition="A",boundaryValue=[0.0d0,0.0d0,500.0d0/1000.0d0*random%gauss(mu=0.0d0,sigma=1.0d0) ])
        t = t + dt
        call siml%solve(dt=dt,timeIntegral="Nemwark-beta",use_same_matrix=.true.)
        
        v = reshape(siml%v,domains(1)%nn(),domains(1)%nd() )
        !if(mod(i_i,10)==0 )then
        !call domains(1)%vtk("v_"+"x_step_"+zfill(i_i/10,5), v(:,1) )
        !call domains(1)%vtk("v_"+"y_step_"+zfill(i_i/10,5), v(:,2) )
        !call domains(1)%vtk("v_"+"z_step_"+zfill(i_i/10,5), v(:,3) )
        !endif
        do j_j=1,37
            call loggers(j_j)%save(t=t)
        enddo
    enddo
    ! destroy
    call siml%remove()
    call mpid%end()
    ![Sim.]
    ! P-wave
    !> 4.6 m/ 35 ms = 131.42 m/s
    ! S-wave
    ! 2.3 m/  35 ms = 65.7 m/s 
    !> 50 m/s
    !Vp/Vs = 4.6/2.3 = 2.0


    ! [Theor]
    !Vp/Vs = sqrt(2.0*(1-0.35)/(1-0.7))
    ! = 2.08

    ! Vp = (101000/(87.30d0/9.80d0)*(1-0.35)/(1.35)/(0.3) )**0.5
    ! = 134.89 m/s
    ! Vs = (101000/(87.300/9.800)/2/1.35)^0.5
    ! =  64.80  m/s

end program main
