Diffuse Scattering from phason for ANC

Calculation of Diffuse Scattering from phason for Al₇₀Ni₁₅Co₁₅ Quasicrystals.

!-------------------------------------------------
! Diffuse Scattering from ANC
!
! Al70Ni15Co15            L=1
!
! ( h1, h2, h3, h4, h5 . h6)
!
!                      by H. Abe
!                      Aug. 26, 1999
!==========================================================
!   a0 = 2.7435
!   as = 2π/a0/sqr(5)    <== a* vector
!
! px(i) = as*COS{2π(i-1)/5} / py(i) = as*SIN{4π(i-1)/5}
!   qx(i) = as*COS{4π(i-1)/5} / qy(i) = as*SIN{4π(i-1)/5}
!
! Ge(1)= h1*px(1)+ h2*px(2) + h3*px(3) + h4*px(4) + h5*px(5)   for external space
! Ge(2)= h1*py(1)+ h2*py(2) + h3*py(3) + h4*py(4) + h5*py(5)
!
! Gi(1)= h1*qx(1)+ h2*qx(2) + h3*qx(3) + h4*qx(4) + h5*qx(5)   for internal space
! Gi(2)= h1*qy(1)+ h2*qy(2) + h3*qy(3) + h4*qy(4) + h5*qy(5)
!
!   G(para) = ( Ge(1), Ge(2) )                    逆格子点
! G(perp) = ( Gi(1), Gi(2) )
!
! Q(para) = G(para) + q(para)       散乱ベクトル
!
!     Sq( , ) : diffuse intensity of calculation
!==========================================================
PROGRAM diffuse
integer, parameter                         :: dp=selected_real_kind(10)
Real,dimension(:,:), allocatable           :: a, b, x, y
Real,dimension(:), allocatable             :: c, fp, fpp
   Real (kind=dp),dimension(:,:), allocatable :: Cinv,Gr,Gl,CGr,Suu,Suv,Svv,sHDS
real (kind=dp),dimension(:), allocatable   :: q, Ge, Gi
real (kind=dp) :: lamda, myu, vvK, uvK, q100, q010, q001, det, S, S1, S2, S3, S4, Ctr
real (kind=dp)      :: THEe, THEi
real                :: dX, dY, WL, g1, g2
integer,dimension(6):: h
integer             :: nAtom, nM, nX, nY
Character (Len=50) :: DR, File_output
Character (len=4)   :: EX
external               MTXINV, Cmatrix, coupling, G_vector
external               Read_ASF, f_open, Fh2_calc
!------------------------------------
! Bragg reflection
!------------------------------------
nM=3         ! number of Matrix for Decagonal phase
ALLOCATE( Cinv(nM*2,nM*2),q(nM),Ge(nM),Gi(nM),Gl(1,nM*2),Gr(nM*2,1), CGr(nM*2,1))
!---------------------------------------------------------
WL=1.493                                   ! Wave Length
! h(1)=3; h(2)=2; h(3)=0; h(4)=0; h(5)=2; h(6)=0 !(6 0 0)    (1 -2 -5 -4.0)phi= 0.0
! h(1)=0; h(2)=2; h(3)=1; h(4)=-1; h(5)=-2; h(6)= 0 !(0 60)    (30 -4-4.0) phi= 90.0
! h(1)=1; h(2)=2; h(3)=2; h(4)=0; h(5)=0; h(6)=0 !(0 3.708 0)(1 0 -3 -3.0) phi= 90.0

h(1)=3; h(2)=2; h(3)=0; h(4)=1; h(5)=3; h(6)=1 !(0.438 4.85 1)(2 0 -4 -3.1)
h(1)=0; h(2)=-1; h(3)=0; h(4)=0; h(5)=0; h(6)=1 !(1.39 0.0 1)(1 0 -1 0.1)
h(1)=-1; h(2)=-2; h(3)=-1; h(4)=0; h(5)=0; h(6)=1 !(3.708 0.0 1)(1 -1 -3 -2.1)
!---------------------------------------------------
call G_vector (h, g1, g2, q100, q010, q001, Ge, Gi, WL)
THEe = ATAN2(Ge(2),Ge(1))
THEi = ATAN2(Gi(2),Gi(1))
!---------------------------------------------------
write(*,'(A16,6I5)')   '   Index     ',h(:)
write(*,'(A16,2F9.4)') '   (hx hy) = ',g1,g2
write(*,'(A16,2F9.4)') 'Ge (a1 a2) = ',Ge(1),Ge(2)
write(*,'(A16,2F9.4)') 'Gi (a3 a4) = ',Gi(1),Gi(2)
write(*,'(A16,2F9.4)') 'θe        = ',THEe
write(*,'(A16,2F9.4)') 'θi        = ',THEi
write(*,'(A16,2F9.4)') '3θe + θi = ',3.0*THEe+THEi

write(*,*) '-------------------------------'
Gl(1,1)=Ge(1); Gl(1,2)=Ge(2) ;Gl(1,3)=Gi(1) ; Gl(1,4)=Gi(2)    ! left
Gr(1,1)=Ge(1); Gr(2,1)=Ge(2) ;Gr(3,1)=Gi(1) ; Gr(4,1)=Gi(2)    ! right
!-------------------------------------
! Read Atomic Scattering Factor
!-------------------------------------
    nAtom = 3     ! 1: Al 2: Ni 3: Co
    ALLOCATE( a(nAtom,4), b(nAtom,4), c(nAtom),fp(nAtom), fpp(nAtom) )
    CALL Read_ASF (a, b, c, fp, fpp, nAtom, WL)
f2ANC=Fh2_calc(Gi,a,b,c,fp,fpp,nAtom)
write(*,'(A16,F9.4)') 'f(cal)^2=',f2ANC
write(*,*) '-------------------'
!-----------------------------------------------------
!   (h, k, 0) plane
!-----------------------------------------------------
nX=31; nY=31
dX=0.02; dY=0.02
! nX=21; nY=21
! dX=0.50; dY=0.50
ALLOCATE( x(nX,nY), y(nX,nY), Suu(nX,nY), Suv(nX,nY), Svv(nX,nY), sHDS(nX,nY) )
CALL Index_hk (g1,g2,X,Y,dX,dY,nX,nY)
!-------------------------------------------------------------
! calculation of the diffuse scattering
!-------------------------------------------------------------
DR='B:\ABE\ANC5\1.493\'
File_output='hk060Ni'
EX='.dat'
File_output=trim(DR)//trim(File_output)//EX
!----------------------------------------------
    lamda = 0.5736     ! Lame constant    [10^12 dyn/cm^2]
myu   = 0.8845    !μ>0, λ+μ>0,μK1>2(K3^2*K4^2)

vvK=0.1      ! phason elastic constant
uvK=0.2         ! coupling constant
!----------------------------------------
write(*,'(A16,F10.5)') ' λ = ',lamda
write(*,'(A16,F10.5)') ' μ = ',myu
write(*,'(A16,F10.5)') 'Kvv = ',vvK
write(*,'(A16,F10.5)') 'Kuv = ',uvK
write(*,*)           File_output
write(*,*) '-------------------'
write(*,*) 'Calculation Start?'
read(*,*) aaaaa
!---------------------------------------------------------
open(8,FILE=file_output)
   write(8,'(2F8.4)') vvK, uvK
   do i=1, nX
     do j=1, nY
    q(1) = (x(i,j)-g1)*q100 * 20.0
    q(2) = (y(i,j)-g2)*q010 * 20.0

    Ctr=0.0
    if (abs(q(1))<1.e-6 .and. abs(q(2))<1.e-6) then
      Suu(i,j)=1.0e6
   Suv(i,j)=0.0
   Svv(i,j)=0.0
   sHDS(i,j)=1.0e6
    else
   sHDS(i,j)=coupling (q, Ge, Gi, lamda, myu, vvK, uvK)
    call Cmatrix (q, lamda, myu, vvK, uvK, Cinv, nM, det)

   if (abs(det)<1.e-6) then
     write(*,*) '-------------------'
           do k = 1, 2*nM
                write(*,'(4F12.6)') Cinv(k,:)
              end do
           write(*,*) '-------------------'
           write(*,'(A6,F12.6)') 'd(M)=',det
     write(*,*) x(i,j), y(i,j)
     read(*,*) aaaa
     Suu(i,j)=1.0e6
   else
     Ctr = Cinv(1,1)+Cinv(2,2)+Cinv(3,3)+Cinv(4,4) ! trace(Cinv)

!------------- total C-matrix) ----------------------
     CGr = matmul(Cinv,Gr)
     S=0.0
           do k = 1, nM*2
             S = S + Gl(1,k)*CGr(k,1)
            end do
!----------------------------------------------------

     S1 = Cinv(1,1)*Ge(1) + Cinv(1,2)*Ge(2) !           (Cuu)^(-1) * G(para)
     S2 = Cinv(2,1)*Ge(1) + Cinv(2,2)*Ge(2) !           (Cuu)^(-1) * G(para)
     Suu(i,j) = Ge(1)*S1 + Ge(2)*S2   ! G(para) * (Cuu)^(-1) * G(para)

     S1 = Cinv(1,3)*Gi(1) + Cinv(1,4)*Gi(2) !           (Cuv)^(-1) * G(perp)
     S2 = Cinv(2,3)*Gi(1) + Cinv(2,4)*Gi(2) !           (Cuv)^(-1) * G(perp)
     S3 = Ge(1)*S1 + Ge(2)*S2           ! G(para) * (Cuu)^(-1) * G(perp)

     S1 = Cinv(3,1)*Ge(1) + Cinv(3,2)*Ge(2) !           (Cuv)^(-1) * G(para)
     S2 = Cinv(4,1)*Ge(1) + Cinv(4,2)*Ge(2) !           (Cuv)^(-1) * G(para)
     S4 = Gi(1)*S1 + Gi(2)*S2           ! G(perp) * (Cuu)^(-1) * G(para)
     Suv(i,j)= S3 + S4

      S1 = Cinv(3,3)*Gi(1) + Cinv(3,4)*Gi(2) !           (Cuv)^(-1) * G(perp)
     S2 = Cinv(4,3)*Gi(1) + Cinv(4,4)*Gi(2) !           (Cuv)^(-1) * G(perp)
     Svv(i,j) = Gi(1)*S1 + Gi(2)*S2   ! G(perp) * (Cuu)^(-1) * G(perp)

end if
end if

    write(8,'(2F9.4, 4F12.1)') q(1), q(2), Suu(i,j)*f2ANC, sHDS(i,j)*f2ANC, Svv(i,j)*f2ANC, S*f2ANC
    write(*,'(2F9.4, 4F12.1)') q(1), q(2), Suu(i,j)*f2ANC, sHDS(i,j)*f2ANC, Svv(i,j)*f2ANC, S*f2ANC
end do
!     read(*,*) aaaaaaa
   end do
close (8)
!----------------------------------
write(*,*) File_output
END PROGRAM diffuse
!****************************************************************************
function coupling (q, Ge, Gi, lamda, myu, vvK, uvK) result(S)
!****************************************************************************
    integer, parameter                      :: dp=selected_real_kind(10)
real (kind=dp),dimension(:),intent(in) :: q(2), Ge(2), Gi(2)
real (kind=dp),intent(in)               :: lamda, myu, vvK, uvK
real (kind=dp)                          :: pi, THE, THEe, THEi, GeAbs, GiAbs, q2
real (kind=dp)                          :: c1, c2, c3, c4, ang1, dTHE, S1, S2
pi = 3.1415926535
THE = ATAN2( q(2), q(1))
THEe = ATAN2(Ge(2),Ge(1))
THEi = ATAN2(Gi(2),Gi(1))

GeAbs = SQRT(Ge(1)**2 + Ge(2)**2)
GiAbs = SQRT(Gi(1)**2 + Gi(2)**2)
q2 = q(1)**2 + q(2)**2
c1=1.0/( vvK-uvK**2/(lamda+2.0*myu) )/q2
c2=1.0/( vvK-uvK**2/ myu )/q2
c3=uvK/( lamda + 2.0*myu )
c4=uvK/myu

ang1= 3.0*THE + THEi
dTHE= THE - THEe

S1 = c1*(GiAbs*cos(ang1)     - c3*GeAbs*cos(dTHE) )**2
   S2 = c2*(GiAbs*sin(ang1)     - c4*GeAbs*sin(dTHE) )**2
! S1 = c1*(GiAbs*cos(3.0*dTHE) - c3*GeAbs*cos(dTHE) )**2   ! 3θe + θi = 0
!   S2 = c2*(GiAbs*sin(3.0*dTHE) - c4*GeAbs*sin(dTHE) )**2   ! 3θe + θi = 0
! S1 = c1*(-GiAbs*cos(3.0*dTHE) - c3*GeAbs*cos(dTHE) )**2   ! 3θe + θi = π
!   S2 = c2*(-GiAbs*sin(3.0*dTHE) - c4*GeAbs*sin(dTHE) )**2   ! 3θe + θi = π
S = S1 + S2
end function coupling
!**********************************************************************************************
SUBROUTINE Index_hk (g1,g2,X,Y,dX,dY,nX,nY)
!**********************************************************************************************
    integer,intent(in)              :: nX, nY
real,intent(in)                 :: dX, dY, g1,g2
real,dimension(:,:),intent(OUT) :: X(nX,nY), Y(nX,nY)
real                            :: Xs, Ys, x1, x2
Xs=g1-(nX-1)/2.0*dX
Ys=g2-(nY-1)/2.0*dY
    x1=Xs; x2=Ys
do i=1,nX
   do j=1,nY
x1=anint(1000*x1)/1000
x2=anint(1000*x2)/1000
X(i,j)=x1; Y(i,j)=x2
!       write(*,*) X(i,j),Y(i,j)
        x1=x1+dX
   end do
   x1=Xs
   x2=x2+dY
!   read(*,*) aaaaa
end do
end subroutine Index_hk
!****************************************************************************
subroutine G_vector (h, g1, g2, q100, q010, q001, Ge, Gi, WL)
!****************************************************************************
    integer, parameter                      :: dp=selected_real_kind(10)
    integer,dimension(:),intent(in)         :: h(6)
    real,intent(in)                         :: WL
    real,intent(out)                        :: g1, g2
real (kind=dp),intent(out)              :: q100, q010
real (kind=dp),dimension(:),intent(out) :: Ge(2), Gi(2)
    real (kind=dp),dimension(:)             :: px(5),py(5),qx(5),qy(5)
    real (kind=dp)                          :: pi, ps, a0, b0, c0, q001,p
    real (kind=dp)                          :: x,y,z,q,ome,a,ar,r,GeABS,GiAbs

    pi = 3.1415926535
    ps = 0.5384385                              ! [ A^(-1) ] p*'=4πsinθ/λ
    p=ps*2.0*cos(pi/10.0)      !             p*= p*'2cos(π/10)
a0 = 8.703785; b0 = 7.391747; c0 = 4.078059 ! [A]
!---------------------------------------------------------------
    q100 = 4.0*pi / a0 / 2.0   ! [ A^(-1) ] q[100]=4πsinθ/λ
    q010 = 4.0*pi / b0 / 2.0
    q001 = 4.0*pi / c0 / 2.0
!---------------------------------------------
    x = 0.0; y = 0.0; z = 0.0
    do i = 1,5
      x = x + h(i) * p * COS(2.0*pi/5.0*(i-1))
      y = y + h(i) * p * SIN(2.0*pi/5.0*(i-1))
    end do
    z = h(6) * q001
Ge(3)=z
q = SQRT(x**2 + y**2 + z**3)                  ! q = 4πsinθ/λ
    ome = asin(q * WL / 4.0 / pi) * 180.0/pi   ! [deg.]
    g1 = x/q100; g2 = y/q010
write(*,'(A16,2F9.4)') '( x, y )= ', g1,g2
    write(*,'(A16,2F9.4)') '(2θ,ω )= ', 2.0*ome, ome
    write(*,'(A16, F9.4)') '4π(sinθ/λ) =', q
!--------------------------------------------
!   p* = 2 p*' cos(π/10)
!   q* = 2 q*' cos(3π/10)
!---------------------------------------------
    a = 2.7435         ! decagonal lattice constant [A]
    ar = 2.0 * pi / a   ! a* of reciprocal space     [A^-1]
    r = ar / SQRT(5.0)
    do i = 1,5
      px(i)=r*COS(2.0*pi*(i-1)/5.0); py(i)=r*SIN(2.0*pi*(i-1)/5.0)
     qx(i)=r*COS(4.0*pi*(i-1)/5.0); qy(i)=r*SIN(4.0*pi*(i-1)/5.0)
end do

Ge(1)=0.0; Ge(2)=0.0; Gi(1)=0.0; Gi(2)=0.0; Gi(3)=0.0
do i=1,5
Ge(1)=Ge(1)+h(i)*px(i); Ge(2)=Ge(2)+h(i)*py(i)
Gi(1)=Gi(1)+h(i)*qx(i); Gi(2)=Gi(2)+h(i)*qy(i)
end do

!   GeAbs = SQRT(Ge(1)**2 + Ge(2)**2)/2.0/COS(    pi/10.0)
!   GiAbs = SQRT(Gi(1)**2 + Gi(2)**2)/2.0/COS(3.0*pi/10.0)
   GeAbs = SQRT(Ge(1)**2 + Ge(2)**2 + Ge(3)**2)
    GiAbs = SQRT(Gi(1)**2 + Gi(2)**2 )

   write(*,'(A16,F9.4)') 'G(external)=',GeAbs
write(*,'(A16,F9.4)') 'G(internal)=',GiAbs
end subroutine G_vector
!****************************************************************************
subroutine Cmatrix (q, lamda, myu, vvK, uvK, Cinv, nM, det)
!****************************************************************************
    integer, parameter                        :: dp=selected_real_kind(10)
real (kind=dp),dimension(:),intent(in)    :: q(nM)
   integer,intent(in)                        :: nM
real (kind=dp),dimension(:,:),intent(out) :: Cinv(nM*2,nM*2)
real (kind=dp),intent(in)                 :: lamda, myu, vvK, uvK
    real (kind=dp),intent(out)                :: det
real (kind=dp),dimension(nM*2,nM*2)       :: Ce,Cc
    real (kind=dp),dimension(nM,nM)           :: Cuu, Cvv, Cuv, Cvu
    real (kind=dp),dimension(nM)    :: uq
real (kind=dp)                  :: q2, qm
real (kind=dp),parameter        :: tau= 1.618033988749895_dp
integer :: n

q2=0.0
do i=1,nM
   q2 = q2+ q(i)**2
end do
qm = sqrt(q2)                  ! magnitude of q on quasiperiodic plane
do i=1,nM
   uq(i)=q(i)/qm      ! unit vectors of q
end do
!------------------------
! Initialize parameters
!------------------------
    do i = 1, nM
   Cuu(i,:)=0._dp; Cvv(i,:)=0._dp; Cuv(i,:)=0._dp; Cvu(i,:)=0._dp
    end do
!--------------------------------------------------------
    Cuu(1, 1) = q2*( (lamda+2.0*myu)*uq(1)**2 + myu*(1-uq(1)**2) )
    Cuu(2, 2) = q2*( (lamda+2.0*myu)*uq(2)**2 + myu*(1-uq(2)**2) )
    Cuu(1, 2) = q2*( (lamda+2.0*myu)*uq(1)*uq(2) - myu*uq(1)*uq(2) )
    Cuu(2, 1) = q2*( (lamda+2.0*myu)*uq(2)*uq(1) - myu*uq(2)*uq(1) )
!--------------------------------------------------------
Cvv(1, 1) = vvK*q2
    Cvv(2, 2) = vvK*q2
Cvv(1, 2) = 0.0_dp
Cvv(2, 1) = 0.0_dp
!--------------------------------------------------------
    Cvu(1,1) = uvK*q2*( uq(1)**2-uq(2)**2 )
    Cvu(2,2) = uvK*q2*( uq(2)**2-uq(1)**2 )
    Cvu(1,2) = -2.0*uvK*q2*uq(1)*uq(2)
    Cvu(2,1) = -2.0*uvK*q2*uq(1)*uq(2)
!-------------------------------------------------------
Cuv=transpose(Cvu)    ! Cuv <== t(Cvu)
!--------------------------------------------------------
    do i = 1, nM
      do j = 1, nM
     Ce(i, j ) = Cuu(i,j)
        Ce(i, j+2) = Cuv(i,j)
     Ce(i+2,j ) = Cvu(i,j)
     Ce(i+2,j+2) = Cvv(i,j)
   end do
end do
!------------------------------------
! write(*,*) '-------------------'
!   do i = 1, nM*2
!      write(*,'(4F12.6)') Ce(i,:)
!    end do
! write(*,*) '-------------------'
! read(*,*) aaaaa
!------------------------------------
   n=nM*2
call MTXINV(Ce,Cinv,n,det)
!------------------------------------
! write(*,*) '-------------------'
!   do i = 1, n
!      write(*,'(4F12.6)') Cinv(i,:)
!    end do
! write(*,*) '-------------------'
! write(*,'(A6,F12.6)') 'd(M)=',det
! read(*,*) aaaaa
!------------------------------------
   Cc=matmul(Ce,Cinv)

! do i = 1, n
!      write(*,'(4F12.6)') Cc(i,:)
!    end do
! write(*,*) '-------------------'
! read(*,*) aaaaa
end subroutine Cmatrix
!************************************************************************
Subroutine MTXINV(a,b,nc,det)
!************************************************************************
    integer,intent(IN)            :: nc
integer,dimension(:)          :: ja(nc),ka(nc)
integer, parameter            :: dp=selected_real_kind(10)
real (kind=dp),dimension(:,:),intent(in) :: a(nc,nc)
real (kind=dp),dimension(:,:),intent(out) :: b(nc,nc)
real (kind=dp) , intent(out) :: det
    real (kind=dp)                :: Amax, work

do i=1,nc
b(i,:)=a(i,:)
end do

det=1._dp
do l=1,nc
   Amax=0._dp
   do j=l,nc
     do k=l,nc
    if ( abs(Amax) <= abs(b(j,k)) ) then
      Amax=b(j,k)
   ja(l)=j; ka(l)=k
    end if
end do
   end do

   if (Amax==0._dp) then
     det=0._dp
     return
   end if

   j=ja(l)
   if (j/=l) then
     do k=1,nc
    work=b(l,k)
    b(l,k)=b(j,k)
    b(j,k)=-work
end do
   end if

   k=ka(l)
   if (k/=l) then
     do j=1,nc
    work=b(j,l)
    b(j,l)=b(j,k)
    b(j,k)=-work
end do
   end if

   do j=1,nc
     if (j/=l) then
    b(j,l)=-b(j,l)/Amax
end if
   end do

   do j=1,nc
     if (j/=l) then
    do k=1,nc
      if (k/=l) then
     b(j,k)=b(j,k)+b(j,l)*b(l,k)
   end if
    end do
end if
   end do

   do k=1,nc
     if (k/=l) then
    b(l,k)=b(l,k)/Amax
end if
   end do

   b(l,l)=1._dp/Amax
   det=det*Amax
end do
!--------------------------
do l=1,nc
   m=nc-l+1
   k=ja(m)
   if(k>m) then
do j=1, nc
    work=b(j,m)
    b(j,m)=-b(j,k)
    b(j,k)=work
end do
   end if

   j=ka(m)
   if(j>m) then
do k=1, nc
    work=b(m,k)
    b(m,k)=-b(j,k)
    b(j,k)=work
end do
   end if
end do
end subroutine MTXINV
!*****************************************************************
Subroutine Read_ASF (a, b, c, fp, fpp, nAtom,WL)
!*****************************************************************
    integer,intent(IN)               :: nAtom
    real,dimension(:,:),intent(out) :: a(nAtom,4), b(nAtom,4)
    real,dimension(:),intent(out)    :: c(nAtom),fp(nAtom),fpp(nAtom)
    real,intent(in)                  :: WL
    integer                          :: nA
    Character (Len=50)               :: DR,File_f0

    if (WL==1.493) then
      fp(1)= 0.19624; fpp(1)=0.231    ! Al
      fp(2)=-5.3956;   fpp(2)=0.4811   ! Ni
      fp(3)=-1.79665; fpp(3)=3.43962 ! Co
    elseif (WL==1.531) then
      fp(1)= 0.20288; fpp(1)=0.24261 ! Al
      fp(2)=-3.245469; fpp(2)=0.5034   ! Ni
      fp(3)=-2.295391; fpp(3)=3.57492 ! Co
    elseif (WL==1.613) then
      fp(1)= 0.21721; fpp(1)=0.26839 ! Al
      fp(2)=-2.29763; fpp(2)=0.55316 ! Ni
      fp(3)=-5.5962;   fpp(3)=0.4756   ! Co
    end if

    DR='B:\QB\X-ray\ASF\'
!-----------------------------------------------------
! Read the Atomic Scattering Factor
    File_f0 =trim(DR)//'Al.f0'
    nA=1
    CALL f_open (File_f0, a, b, c, nAtom, nA)

    File_f0 =trim(DR)//'Ni.f0'
    nA=2
    CALL f_open (File_f0, a, b, c, nAtom, nA)

    File_f0 =trim(DR)//'Co.f0'
    nA=3
    CALL f_open (File_f0, a, b, c, nAtom, nA)

!   do i=1,nA
!     do j=1,4
!       write(*,*) a(i,j),b(i,j)
!   end do
!     write(*,*) c(i)
!   write(*,*) fp(i),fpp(i)
!   read(*,*) aaaaa
!   end do
end Subroutine Read_ASF
!********************************************************************************
FUNCTION Fh2_calc(Gi,a,b,c,fp,fpp,nAtom) result(Fh2)
!********************************************************************************
    integer, parameter                     :: dp=selected_real_kind(10)
real (kind=dp),dimension(:),intent(IN) :: Gi(2)
integer,intent(IN)              :: nAtom
real,dimension(:,:),intent(IN) :: a(nAtom,4), b(nAtom,4)
real,dimension(:),intent(IN)    :: c(nAtom),fp(nAtom),fpp(nAtom)
    real,dimension(10)              :: Fr
Real                            :: Freal, Fimag
real (kind=dp)                  :: pi,g
pi=3.1415926535_dp

g = SQRT(Gi(1)**2 + Gi(2)**2)/2.0/COS(3.0*pi/10.0)/2/pi   ! g = sinθ/λ
!---------------------------
!   Atomic Scattering Factor
!--------------------------------
!     1: Al / 2: Ni / 3: Co
!--------------------------------
do k=1,nAtom
   Fr(k)=0.0
   do l=1,4
Fr(k) = Fr(k) + a(k,l)*exp( -b(k,l)*g**2 )
   end do
   Fr(k) = Fr(k) + c(k) + fp(k)
end do
!   average Structure Factor for ANC
Freal=0.70*Fr(1) + 0.15*Fr(2) + 0.15*Fr(3)
Fimag=0.70*fpp(1)+ 0.15*fpp(2) +0.15*fpp(3)
    write(*,'(A16,2F9.4)') 'F(real imag)=',Freal,Fimag

Fh2 = Freal**2 + Fimag**2
end function Fh2_calc
!************************************************************************
subroutine f_open (file_f0, a, b, c, nAtom, nA)
!************************************************************************
integer,intent(in)              :: nAtom, nA
    character (Len=50) ,intent(in) :: File_f0
real,dimension(:,:),intent(OUT) :: a(nAtom,4), b(nAtom,4)
real,dimension(:),intent(OUT)   :: c(nAtom)

open( 2,FILE=file_f0 )
   do i=1,4
read(2,*) a(nA,i)
read(2,*) b(nA,i)
   end do
     read(2,*) c(nA)
close(2)
end subroutine f_open
!************************************************************************

References

H. Abe et al., to be published.

Back to ABE

Back to Department of Materials Science and Engineering

Back to research

Department of Materials Science and Engineering
National Defense Academy

Last Modified: April 1, 2009

Calculation of Diffuse Scattering from phason for Al70Ni15Co15 Quasicrystals.

Calculation of Diffuse Scattering from phason for Al₇₀Ni₁₅Co₁₅ Quasicrystals.