[015]結晶粒の境界

 投稿者：mike  投稿日：2017年 5月 1日(月)09時33分22秒
  Morseポテンシャルによる結晶粒の境界（粒界）を観察するための分子動力学法プログラムを公開します。
結晶には方位があり、方位のそろわない境界には不整合が存在します。この不整合は材料の物性に影響するため、
分子動力学法による研究が盛んな分野の一つです。
本プログラムでは、２つの方位の異なる結晶をくっつけて粒界をつくります。
粒界は不安定で、時間とともに動きます。

表示の説明：
結合の色は方位を表します。60度ごとに同じ色になるようにすることにより、
1つの結晶が同じ色に見えるようにしています。

試験環境：
本プログラムは十進BASIC 6.6.2.2 / macOS 10.7.5, 十進BASIC Ver 7.8.0 / windows 10　でテストしました。

-------------

!
! ========= molecular dynamics 2D - Morse potential ==========
!
! 015grainBoundaryMMD2D.bas
!   Copyright(C) 2017 Mitsuru Ikeuchi
!
! ver 0.0.1   2017.04.18  created
!
OPTION ARITHMETIC NATIVE
DECLARE EXTERNAL SUB mmd2d.setInitialCondition, mmd2d.moveParticles, mmd2d.drawParticles
DECLARE EXTERNAL FUNCTION INKEY$
LET tempMode = 1    !tempMode: 0:adiabatic  1:constant temperature
LET contTemp = 300  !contTemp: controled constant temperature(K)
LET drawMode = 1    !drawMode: 0:bond  1:circle+bond 2:velocitySpace
DIM menu$(8)
LET menu$(1) = "Fe 2.8 6x6"
LET menu$(2) = "Fe 4.0 6x6"
LET menu$(3) = "Cu 4.0 6x6"
LET menu$(4) = "Al 4.0 6x6"
LET menu$(5) = "Ca 4.0 6x6"
LET menu$(6) = "Ba 4.0 6x6"
LET menu$(7) = "W  4.0 6x6"
LET menu$(8) = "continue"

!setInitialCondition(molKind,boxSizeInNM,xtalSizeInNM,contTemp)
CALL setInitialCondition(3,6.0,4.0,contTemp) !molKind: 3:Fe  5:Al  13:Ba  15:K  18:Xe

DO
   CALL moveParticles(tempMode,contTemp)
   CALL drawParticles(tempMode,contTemp,drawMode)
   LET S$=INKEY$
   IF S$="." THEN
      EXIT DO
   ELSEIF S$="0" THEN
      LOCATE VALUE (1),RANGE 0 TO 4000 : temp
      IF temp>1 THEN
         LET tempMode = 1
         LET contTemp = temp
      ELSE
         LET tempMode = 0
      END IF
   ELSEIF S$="1" THEN
      LET drawMode = MOD(drawMode+1,3)
   ELSEIF S$="9" THEN
      LOCATE CHOICE (menu$) :nmenu
      IF nmenu=1 THEN      !--- Fe 2.8 box=6x6nm
         CALL setInitialCondition(3,6.0,2.8,contTemp)
      ELSEIF nmenu=2 THEN  !--- Fe 4.0
         CALL setInitialCondition(3,6.0,4.0,contTemp)
      ELSEIF nmenu=3 THEN  !--- Cu 4.0
         CALL setInitialCondition(7,6.0,4.0,contTemp)
      ELSEIF nmenu=4 THEN  !--- Al 4.0
         CALL setInitialCondition(5,6.0,4.0,contTemp)
      ELSEIF nmenu=5 THEN  !--- Ca 4.0
         CALL setInitialCondition(9,6.0,4.0,contTemp)
      ELSEIF nmenu=6 THEN  !--- Ba 4.0
         CALL setInitialCondition(11,6.0,4.0,contTemp)
      ELSEIF nmenu=7 THEN  !--- W 4.0
         CALL setInitialCondition(0,6.0,4.0,contTemp)
      ELSEIF nmenu=8 THEN  !contine
      !
      END IF

   END IF
LOOP

END

EXTERNAL FUNCTION INKEY$ !from decimal BASIC library
OPTION ARITHMETIC NATIVE
SET ECHO "OFF"
LET S$=""
CHARACTER INPUT NOWAIT: S$
LET INKEY$=S$
END FUNCTION

! ---------- molecular dynamics 2D - Morse potential ----------
!
!  method: velocity Verlet algorithm
!    (1) vi = vi + (Fi/mi)*(0.5dt)
!    (2) ri = ri + vi*dt
!    (3) calculation Fi <- {r1,r2,...,rn} Fi=sum(Fij,j=1 to n),Fij=F(ri-rj)
!    (4) vi = vi + (Fi/mi)*(0.5dt)
!    (6) goto (1)
!  potential: Morse V(r) = D*((1-EXP(-A*(r-r0)))^2-1)
!                        = D*(EXP(-2*A*(r-r0))-2*EXP(-A*(r-r0)))
!    (D:dissociation energy, r0:bond length, A:width parameter { A=SQR(k/(2*D)) }
!             force F(r) = -dV(r)/dr
!                        = 2*D*A*y*(y-1), y=EXP(-A*(r-r0))
MODULE mmd2d
MODULE OPTION ARITHMETIC NATIVE
PUBLIC SUB setInitialCondition !(molKind,boxSizeInNM,xtalSizeInNM,contTemp)
PUBLIC SUB moveParticles !(tempMode,contTemp)
PUBLIC SUB drawParticles !(drawMode:0:realSpace 1:velocitySpace)
SHARE NUMERIC molecKind, sysTime, dt, nMolec, xMax, yMax, mass, Dmrs, Amrs, r0mrs, rCutoff, hh
SHARE NUMERIC xx(500),yy(500)           ! (xx(i),yy(i))  : position of i-th particle
SHARE NUMERIC vx(500),vy(500)           ! (vx(i),vy(i))  : velocity of i-th particle
SHARE NUMERIC ffx(500),ffy(500)         ! (ffx(i),ffy(i)): total force of i-th particle
SHARE NUMERIC reg(500,0 TO 100)         ! register near molec reg(i,0):number of near i-th molec
SHARE NUMERIC molecData(0 TO 20,0 TO 3) ! molecule 0:mass, 1:epsilon, 2:sigma, 3:dt
SHARE STRING molecSTR$(0 TO 20)         ! molec string
SHARE NUMERIC forceTable(0 TO 1200)     ! force table
LET molecKind = 3          ! 3:Fe
LET sysTime = 0.0          ! system time (s) in the module
LET dt = 5.0*1.0e-15       ! time step (s)
LET nMolec = 100           ! number of particles
LET xMax = 6.0E-9          ! x-Box size (m)
LET yMax = 6.0E-9          ! y-Box size (m)
LET mass = 55.847*1.661e-27! mass of Fe (kg)
LET Dmrs = 0.4174*1.602e-19! D of Morse potential (J)   : energy of dissociation
LET Amrs = 1.3885e10       ! A of Morse potential (1/m) : width parameter
LET r0mrs = 2.845e-10      ! r0 of Morse potential (m)  : bond length
LET rCutoff = 1.0e-9       ! force cutoff radius (m)
LET hh = 1e-12             ! force table step (m)
LET molecData(0,0)=0       ! if molecData(0,0)=0 then set molecData(,)

! ---------- set initial condition

EXTERNAL SUB setInitialCondition(molKind,boxSizeInNM,xtalSizeInNM,contTemp)
   DECLARE EXTERNAL SUB setMolecData,setMolecules,ajustVelocity,setForceTable
   DECLARE EXTERNAL FUNCTION setCrystalBlock
   RANDOMIZE
   CALL setMolecData
   LET molecKind = molKind
   LET mass = molecData(molecKind,0)
   LET Dmrs = molecData(molecKind,1)
   LET Amrs = molecData(molecKind,2)
   LET r0mrs = molecData(molecKind,3)
   LET sysTime = 0.0
   LET xMax = boxSizeInNM*1.0e-9
   LET yMax = boxSizeInNM*1.0e-9
   ! set particles
   LET s = 0.5*(boxSizeInNM-xtalSizeInNM)*1.0e-9
   LET a = xtalSizeInNM*1.0e-9
   LET d = 0.5*r0mrs
   LET ii= setCrystalBlock(1, s, s, 0.5*a-0.5*d, a, PI/2)
   LET nMolec = setCrystalBlock(ii+1, s+0.5*a+0.5*d, s, 0.5*a-0.5*d, a, 0)
   CALL ajustVelocity(contTemp)
   LET rCutoff = MIN(1.0e-9, 3.0*r0mrs)
   CALL setForceTable
   ! set window
   SET WINDOW 0,500,0,500
END SUB

EXTERNAL SUB setMolecData
! Morse potential data
!       0:mass(AU) 1:D(eV)  2:A(m^-1)  3:r0(m)
   DATA   183.85 , 0.9906, 1.4116e10, 3.032e-10 !  0 W
   DATA    95.94 , 0.8032, 1.5079e10, 2.976e-10 !  1 Mo
   DATA    51.996, 0.4414, 1.5721e10, 2.754e-10 !  2 Cr
   DATA    55.847, 0.4174, 1.3885e10, 2.845e-10 !  3 Fe
   DATA    58.71 , 0.4205, 1.4199e10, 2.780e-10 !  4 Ni
   DATA    26.98 , 0.2703, 1.1646e10, 3.253e-10 !  5 Al
   DATA   207.19 , 0.2348, 1.1836e10, 3.733e-10 !  6 Pb
   DATA    63.54 , 0.3429, 1.3588e10, 2.866e-10 !  7 Cu
   DATA   107.87 , 0.3323, 1.3690e10, 3.115e-10 !  8 Ag
   DATA    40.08 , 0.1623, 0.8054e10, 4.569e-10 !  9 Ca
   DATA    87.62 , 0.1513, 0.7878e10, 4.988e-10 ! 10 Sr
   DATA   137.34 , 0.1416, 0.6570e10, 5.373e-10 ! 11 Ba
   DATA    22.99 , 0.0633, 0.5900e10, 5.336e-10 ! 12 Na
   DATA    39.102, 0.0542, 0.4977e10, 6.369e-10 ! 13 K
   DATA    85.47 , 0.0464, 0.4298e10, 7.207e-10 ! 14 Rb
   DATA   132.905, 0.0449, 0.4157e10, 7.557e-10 ! 15 Cs
   DATA    20.183, 0.0031, 1.6500e10, 3.076e-10 ! 16 Ne
   DATA    39.948, 0.0104, 1.3400e10, 3.816e-10 ! 17 Ar
   DATA    83.80 , 0.0141, 1.2500e10, 4.097e-10 ! 18 Kr
   DATA   131.30 , 0.0200, 1.2400e10, 4.467e-10 ! 19 Xe
   DATA   200.59 , 0.0734, 1.4900e10, 3.255e-10 ! 20 Hg

   DATA "W" ,"Mo","Cr","Fe","Ni","Al","Pb","Cu","Ag","Ca","Sr"
   DATA "Ba","Na","K" ,"Rb","Cs","Ne","Ar","Kr","Xe","Hg"

   IF molecData(0,0)=0 THEN
      MAT READ molecData
      FOR i=0 TO 20
         LET molecData(i,0) = molecData(i,0)*1.661e-27 !mass(AU)--> (kg)
         LET molecData(i,1) = molecData(i,1)*1.602e-19 !eps(eV) --> (J)
      NEXT i
      MAT READ molecSTR$
   END IF
END SUB

EXTERNAL function setCrystalBlock(ii, x0, y0, xLen, yLen, theta)
   DECLARE EXTERNAL SUB setParticle
   LET iip = ii
   LET a = 0.98*r0mrs
   LET b = 0.866025*a
   LET leng = xLen
   IF (leng<yLen) THEN LET leng = yLen
   LET leng = 1.5*leng
   LET nx = INT(leng/b) + 1
   LET ny = INT(leng/a) + 1
   LET sth = SIN(theta)
   LET cth = COS(theta)
   FOR i=0 TO nx-1
      LET x = b*i - leng/2.0
      FOR j=0 TO ny-1
         LET y = a*j - leng/2.0
         IF MOD(i,2)=1 THEN LET y = y + 0.5*a
         LET xp = x0 + xLen/2.0 + cth*x - sth*y
         LET yp = y0 + yLen/2.0 + sth*x + cth*y
         IF (xp>=x0 AND xp<=x0+xLen AND yp>=y0 AND yp<=y0+yLen) THEN
            CALL setParticle(iip, xp, yp)
            LET iip = iip + 1
         END IF
      NEXT j
   NEXT i
   LET setCrystalBlock = iip - 1
end function

EXTERNAL SUB setParticle(i, x, y)
   LET xx(i) = x
   LET yy(i) = y
   LET vx(i) = 200.0*(RND+RND+RND+RND+RND+RND-3)
   LET vy(i) = 200.0*(RND+RND+RND+RND+RND+RND-3)
   LET ffx(i) = 0.0
   LET ffy(i) = 0.0
END SUB

EXTERNAL SUB setForceTable
   DECLARE EXTERNAL  FUNCTION cutoff
   FOR ir=10 TO 1200
      LET r = ir*hh
      LET y = EXP(-Amrs*(r-r0mrs))
      LET forceTable(ir) = cutoff(r)*2.0*Dmrs*Amrs*y*(y-1)
   NEXT ir
   FOR ir=0 TO 9
      LET forceTable(ir) = forceTable(10)
   NEXT ir
END SUB

EXTERNAL FUNCTION cutoff(r)
   IF r>0 AND r<0.8*rCutoff THEN
      LET ret = 1
   ELSEIF r>=0.8*rCutoff AND r<rCutoff THEN
      LET ret = 0.5+0.5*COS(PI*(r-0.8*rCutoff)/(0.2*rCutoff))
   else
      LET ret = 0
   END IF
   LET cutoff = ret
END FUNCTION

! ---------- move particles

EXTERNAL SUB moveParticles(tempMode,contTemp) !tempMode 0:adiabatic  1:constant-temp
   DECLARE EXTERNAL SUB moveParticlesDT,ajustVelocity,registerNearMolec
   IF (tempMode=1) THEN CALL ajustVelocity(contTemp)
   CALL registerNearMolec
   FOR i=1 TO 20
      CALL moveParticlesDT
   NEXT i
END SUB

EXTERNAL SUB moveParticlesDT ! velocity Verlet method
   DECLARE EXTERNAL SUB calcForce
   LET a = 0.5*dt/mass
   FOR i=1 TO nMolec
   !LET a = 0.5*dt/mass(i)
      LET vx(i) = vx(i)+a*ffx(i)
      LET vy(i) = vy(i)+a*ffy(i)
      LET xx(i) = xx(i)+vx(i)*dt
      LET yy(i) = yy(i)+vy(i)*dt
   NEXT i
   CALL calcForce
   FOR i=1 TO nMolec
   !LET a = 0.5*dt/mass(i)
      LET vx(i) = vx(i)+a*ffx(i)
      LET vy(i) = vy(i)+a*ffy(i)
   NEXT i
   LET sysTime=sysTime+dt
END SUB

EXTERNAL SUB calcForce
   DECLARE EXTERNAL FUNCTION force,boundaryForce
   LET s = 0.5*3.418e-10
   FOR i=1 TO nMolec
      LET ffx(i) = 0.0
      LET ffy(i) = 0.0
   NEXT i
   FOR i=1 TO nMolec-1
      FOR k=1 TO reg(i,0)-1
         LET j = reg(i,k)
         LET xij = xx(i)-xx(j)
         LET yij = yy(i)-yy(j)
         LET rij = SQR(xij*xij+yij*yij)
         IF rij<rCutoff THEN
            LET f = force(rij)
            LET fxij = f*xij/rij
            LET fyij = f*yij/rij
            LET ffx(i) = ffx(i)+fxij
            LET ffy(i) = ffy(i)+fyij
            LET ffx(j) = ffx(j)-fxij
            LET ffy(j) = ffy(j)-fyij
         END IF
      NEXT k
   NEXT i
   FOR i=1 TO nMolec  ! boundary force
      LET ffx(i) = ffx(i)+boundaryForce(xx(i)+s)+boundaryForce(xx(i)-xMax-s)
      LET ffy(i) = ffy(i)+boundaryForce(yy(i)+s)+boundaryForce(yy(i)-yMax-s)
   NEXT i
END SUB

EXTERNAL FUNCTION force(r) !force(r) <-- forceTable - linear interporation
   LET ir = INT(r/hh)
   LET a = r - ir*hh
   LET force = ((hh-a)*forceTable(ir) + a*forceTable(ir+1))/hh
END FUNCTION

EXTERNAL FUNCTION boundaryForce(r)
   LET r6 = (3.418e-10/r)^6
   LET boundaryForce = (24.0*(0.5*1.711e-21)*r6*(2.0*r6-1.0)/r)
END FUNCTION

EXTERNAL SUB registerNearMolec
   LET rCut = rCutoff+20*2000*dt
   LET rcut2 = rCut*rCut
   FOR i=1 TO nMolec-1
      LET k = 1
      FOR j=i+1 TO nMolec
         LET r2 = (xx(i)-xx(j))*(xx(i)-xx(j))+(yy(i)-yy(j))*(yy(i)-yy(j))
         IF (r2<rcut2) THEN
            LET reg(i,k) = j
            LET k = k + 1
         END IF
      NEXT j
      LET reg(i,0) = k
   NEXT i
END SUB

EXTERNAL FUNCTION maxNearMolec
   LET mx = 0
   FOR i=1 TO nMolec-1
      IF mx<reg(i,0) THEN LET mx = reg(i,0)
   NEXT i
   LET maxNearMolec = mx-1
END FUNCTION

! ---------- utility

EXTERNAL FUNCTION systemTemprature
   LET kB = 1.38e-23 ! Boltzman's constant (J/K)
   LET ek= 0.0       !kinetic energy (J)
   FOR i=1 TO nMolec
      LET ek = ek + 0.5*mass*(vx(i)^2+vy(i)^2)
   NEXT i
   LET systemTemprature = ek/(nMolec*kB) !2D: E/N=kT, 3D: E/N=(3/2)kT
END FUNCTION

EXTERNAL SUB ajustVelocity(temp)
   DECLARE EXTERNAL FUNCTION systemTemprature
   LET r = sqr(temp/systemTemprature)
   FOR i=1 TO nMolec
      LET vx(i) = r*vx(i)
      LET vy(i) = r*vy(i)
   NEXT i
END SUB

! ---------- draw particles

EXTERNAL SUB drawParticles(tempMode,contTemp,drawMode)
   DECLARE EXTERNAL FUNCTION systemTemprature,maxNearMolec
   DECLARE EXTERNAL sub realSpace,velocitySpace,plotBond,plotBondDirection
   SET DRAW MODE HIDDEN
   CLEAR
   IF drawMode=0 OR drawMode=1 THEN !--- 0:circle  1:circle+bond
      call plotBondDirection(drawMode) !plotBond(drawMode)
   ELSEIF drawMode=2 THEN
      call velocitySpace
   END IF

   !--- draw caption
   SET TEXT HEIGHT 10
   SET TEXT COLOR 1 ! black
   LET tmp$ = "adiabatic   constantTemp  "
   PLOT TEXT, AT  50, 70 ,USING "time =#####.## (ps)   temp =####.## (K)":sysTime*1E12,systemTemprature
   PLOT TEXT, AT  50, 55 ,USING molecSTR$(molecKind)&"  N =####":nMolec
   PLOT TEXT, AT  50, 40 ,USING "tempMode = ## ":tempMode
   PLOT TEXT, AT 200, 40 :tmp$(tempMode*12+1:tempMode*12+12)
   PLOT TEXT, AT  50, 25 ,USING "controled Temperature =####.# (K)":contTemp
   PLOT TEXT, AT  50, 10 :"2-dimensional molecular dynamics"
   PLOT TEXT, AT  50,480 :"[.] exit      [0] temp > value "
   PLOT TEXT, AT  50,460 :"[1] dispMode  [9] menu > select"
   SET DRAW MODE EXPLICIT
END SUB

EXTERNAL sub plotBond(drawMode)
   LET boxSize = 300
   LET mag = boxSize/xMax
   LET xp = 100
   LET yp = 100
   SET LINE COLOR 1 ! black : !--- box
   PLOT LINES: xp,yp; boxSize+xp,yp; boxSize+xp,boxSize+yp; xp,boxSize+yp; xp,yp
   SET TEXT HEIGHT 6
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT xp,boxSize+2+yp ,USING  "box size =##.# x ##.# (nm)":xMax*1e9,yMax*1e9
   FOR i=1 TO nMolec
      SET LINE COLOR 8 ! gray
      DRAW circle WITH SCALE(r0mrs/2.0*mag)*SHIFT(xx(i)*mag+xp,yy(i)*mag+yp)
   NEXT i
   IF drawMode=1 THEN
      FOR i=1 TO nMolec-1
         FOR k=1 TO reg(i,0)-1
            LET j = reg(i,k)
            LET r = SQR((xx(i)-xx(j))*(xx(i)-xx(j))+(yy(i)-yy(j))*(yy(i)-yy(j)))
            IF r<1.2*r0mrs THEN
               IF r<0.93*r0mrs THEN
                  SET LINE COLOR 4 !red 13 !'oleave
               ELSEIF r<1.03*r0mrs THEN
                  SET LINE COLOR 3 !green
               ELSEIF r<1.08*r0mrs THEN
                  SET LINE COLOR 2 !blue
               ELSE
                  SET LINE COLOR 8 !gray
               END IF
               PLOT LINES: xx(i)*mag+xp,yy(i)*mag+yp;xx(j)*mag+xp,yy(j)*mag+yp
            END IF
         NEXT k
      NEXT i
   END IF
END sub

EXTERNAL SUB plotBondDirection(drawMode)
   LET boxSize = 300
   LET mag = boxSize/xMax
   LET xp = 100
   LET yp = 100
   SET COLOR MIX(240) 1.0,0.2,0.2
   SET COLOR MIX(241) 0.6,0.6,0.0
   SET COLOR MIX(242) 0.0,1.0,0.0
   SET COLOR MIX(243) 0.0,0.6,0.6
   SET COLOR MIX(244) 0.2,0.2,1.0
   SET COLOR MIX(245) 0.8,0.0,0.8
   SET LINE COLOR 1 ! black : !--- box
   PLOT LINES: xp,yp; boxSize+xp,yp; boxSize+xp,boxSize+yp; xp,boxSize+yp; xp,yp
   SET TEXT HEIGHT 6
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT xp,boxSize+2+yp ,USING  "box size =##.# x ##.# (nm)":xMax*1e9,yMax*1e9
   FOR i=1 TO nMolec
      SET LINE COLOR 8 ! gray
      DRAW circle WITH SCALE(r0mrs/2.0*mag)*SHIFT(xx(i)*mag+xp,yy(i)*mag+yp)
   NEXT i
   IF drawMode=1 THEN
      FOR i=1 TO nMolec-1
         FOR j=i TO nMolec
            LET r = SQR((xx(i)-xx(j))*(xx(i)-xx(j))+(yy(i)-yy(j))*(yy(i)-yy(j)))
            IF r<1.2*r0mrs THEN
               LET yij = yy(i)-yy(j)
               IF yij=0.0 THEN LET yij = 1e-20
               LET th = 3.0*(ATN((xx(i)-xx(j))/yij)+PI/2.0)/PI
               LET col = 240 + INT((th-INT(th))*6)
               SET LINE COLOR col
               PLOT LINES: xx(i)*mag+xp,yy(i)*mag+yp;xx(j)*mag+xp,yy(j)*mag+yp
            END IF
         NEXT j
      NEXT i
   END IF
END sub

EXTERNAL sub velocitySpace
   LET boxSize = 300
   LET xp = 100
   LET yp = 100
   SET LINE COLOR 1 !black : axis
   PLOT LINES: xp,boxSize/2+yp; boxSize+xp,boxSize/2+yp !vx-axis
   PLOT LINES: boxSize/2+xp,yp; boxSize/2+xp,boxSize+yp !vy-axis
   SET TEXT HEIGHT 6
   SET TEXT COLOR 1 ! black
   PLOT TEXT, AT boxSize+xp,boxSize/2+yp: "vx"
   PLOT TEXT, AT boxSize+xp,boxSize/2-12+yp: "2000m/s"
   PLOT TEXT, AT boxSize/2-12+xp,boxSize+yp: "vy  2000m/s"
   PLOT TEXT, AT boxSize/2-8+xp,boxSize/2-10+yp: "0"
   PLOT TEXT, AT xp,boxSize+8+yp: "velocity space (vx,vy)"
   LET mag = boxSize/4000
   FOR i=1 TO nMolec
      SET LINE COLOR 2 ! blue
      DRAW circle WITH SCALE(5)*SHIFT(vx(i)*mag+boxSize/2+xp,vy(i)*mag+boxSize/2+yp)
   NEXT i
END sub

END MODULE
 

戻る