Devil Maker: Tokyo_はじめに

たまには, 息抜きにゲームでも.

日本のソシャゲは, どれも最終的にはガチャで良いものを引かないと強くなれない飽和状態によくなる.
また, ガチャを引くアイテムも, ゲームが進行するごとに手に入り難く, 課金を遠回しに要求している.

お金を使わないと強くなれないというのは, ゲームとしてどうなのだろうか.

今回は, 以前からプレイしているゲームを紹介する.

DevilMaker: Tokyoは, 神話などをモチーフにしたカードゲームである.
海外のソーシャルゲームであり. 韓国版と英語版のものがある.
日本語版も以前はあったが, 現在は運営されてない.

おすすめする理由

1, ガチャが重要視されてない.
このゲームは, ガチャで☆3~5のカードがドロップするが, 最高レア度は☆7である.
課金したからといって, 最高レア度を入手することはできず,強化して行けばどれも強いカードになる.
また, ガチャを引くためのアイテムも, 購入せずともゲーム内でどんどん手に入るので困ることは無い.
(リセマラは重要でない)

2, 課金は主にプレイ時間の短縮を目的としている.
課金の主な内容が, ゲーム内アイテムの販売である. これは普通にプレイしていれば入手できるものばかりである.
他にも経験値増加やドロップ率増加アイテムも販売している. これらも, ランダムでダンジョンで同じ効果が出現する.
アイテム収集する時間, ゲームをずっとプレイする時間はないけど, ゲームは楽しみたい人用に, 課金アイテムが存在する.
つまり, このゲームは, 課金はプレイ時間を買うことに相当する.
たまに有料のイベントガチャもあるが, 最強カードがドロップなんてことは, 無い.
(イベントで上位ランクに入れば, ゲーム内課金通貨も報酬で手に入る. )

3, 好きなカードがほぼ確実に手に入る.
たとえランクの低いカードを手に入れても強化・進化をして行けば, どれも最高ランクまで強化可能である.
また, 毎月好きな☆5カードが手に入るシステムもあるので, 物欲センサーがあろうとなかろうと関係ない.

4, 絵が綺麗
正直, 適当な絵だったら, プレイしていない.

ゲームが外国語というのはやや抵抗があるだろうが, ゲーム自体は難しくないので(艦これみたいに半自動進行)すぐになれるだろう.

これから, このゲームについて, 少しずつだが, 説明していく.

二重振り子

座標の取り方は下図のように取ります。棒の伸び縮みは無いものとします。
二重振り子座標
どういう解き方でもいいですが、ここでは

  1. デカルト座標\(L(x,y)\)でラグランジアンを記述
  2. デカルト座標から座標変換し、\((r,\theta)\)でラグランジアンを記述
  3. 新たな座標系で運動方程式を導く

という順で解いていきます。

[adsense1]

1, デカルト座標でのラグランジアンLは(運動エネルギーK)-(位置エネルギーU)と書けるため、
\(
L(x_1,\dot{x}_1,y_1,\dot{y}_1,x_2,\dot{x}_2,y_2,\dot{y}_2)\\
\displaystyle =\frac{1}{2}m_1(\dot{x}_1^2+\dot{y}_1^2)+\frac{1}{2}m_1(\dot{x}_1^2+\dot{y}_2^2)-(-m_1gy_1-m_2gy_2)
\)
と書けます。

2, デカルト座標から座標変換
式を簡単にするために座標変換を行います。新しい座標\((r_1,\theta_1,r_2,\theta_2)\)とデカルト座標\((x_1,y_1,x_2,y_2)\)の関係式は
\(
\begin{align}
x_1&=r_1\sin{\theta_1}\\
y_1&=-r_1\cos{\theta_1}\\
x_2&=r_1\sin{\theta_1}+r_2\sin{\theta_2}\\
y_2&=-r_1\cos{\theta_1}-r_2\cos{\theta_2}
\end{align}
\)
という関係があります。各々を時間で微分すれば、
\(
\begin{align}
\dot{x}_1&=\dot{r}_1\sin{\theta_1}+r_1\dot{\theta}_1\cos{\theta_1}\\
\dot{y}_1&=-\dot{r}_1\cos{\theta_1}+r_1\dot{\theta}_1\sin{\theta_1}\\
\dot{x}_2&=\dot{r}_1\sin{\theta_1}+r_1\dot{\theta}_1\cos{\theta_1}+\dot{r}_2\sin{\theta_2}+r_2\dot{\theta}_2\cos{\theta_2}\\
\dot{y}_2&=-\dot{r}_1\cos{\theta_1}+r_1\dot{\theta}_1\sin{\theta_1}-\dot{r}_2\cos{\theta_2}+r_2\dot{\theta}_2\sin{\theta_2}
\end{align}
\)
これらをラグランジアン\(L(x_1,\dot{x}_1,y_1,\dot{y}_1,x_2,\dot{x}_2,y_2,\dot{y}_2)\)に代入します。すると、新たな座標系でのラグランジアン\(L(r_1,\dot{r}_1,\theta_1,\dot{\theta}_1,r_2,\dot{r}_2,\theta_2,\dot{\theta}_2)\)が得られます。
\(
\begin{align}
L(r_1,\dot{r}_1,\theta_1,\dot{\theta}_1,& r_2,\dot{r}_2,\theta_2,\dot{\theta}_2) \\
=&\frac{1}{2}m_1(\dot{r}_1^2+r_1^2\dot{\theta}_1^2)+\frac{1}{2}m_2\left[\dot{r}_1^2+\dot{r}_2^2+r_1^2\dot{\theta}_1^2+r_2^2\dot{\theta}_2^2 \right. \\
&\left.+2(\dot{r}_1 r_2 \dot{\theta}_2-r_1\dot{r}_2 \dot{\theta}_1)\sin{(\theta_1-\theta_2)}+2(\dot{r}_1 \dot{r}_2 +r_1 r_2 \dot{\theta}_1 \dot{\theta}_2)+\cos{(\theta_1-\theta_2)}\right] \\
&+m_1gr_1\cos{\theta_1}+m_2g(r_1\cos{\theta_1}+r_2\cos{\theta_2})
\end{align}
\)

僕は先ほど式を簡単にするために座標変換をする、といいました。しかし、新しい座標系におけるラグランジアンはどう見ても元のデカルト座標系のラグランジアンに比べて複雑です。この理由は物理的な意味から来ています。
振り子をつないでいる棒が伸び縮みしないとすると系の自由度は角度\(\theta_1,\theta_2\)の2つです。
となると運動方程式は最高で2本の独立した方程式になるはずです。
しかし、デカルト座標の場合うまく自由度を落とすことができず、運動方程式は4つになってしまいます。
そこで棒が伸び縮みを簡単に表すことができる座標系に移ることでうまく方程式の数を減らせます。

新しい座標系でのラグランジアンで棒の伸び縮みがないという条件を表すには
\(
\begin{align}
r_1&=l_1\ (l_1\mbox{は定数}) \\
r_2&=l_2\ (l_2\mbox{は定数})
\end{align}
\)
と書けるわけで、また、
\(
\begin{align}
\dot{r}_1&=0 \\
\dot{r}_2&=0
\end{align}
\)
となるわけです。

\(m_1=m_2=m,\ l_1=l_2=l\)という場合を特に考えると、ラグランジアンは
\(
\displaystyle L(\theta_1,\dot{\theta}_1,\theta_2,\dot{\theta}_2)=ml^2\left[\dot{\theta}_1^2+\frac{1}{2}\dot{\theta}_2^2+\dot{\theta}_1\dot{\theta}_2\cos{(\theta_1-\theta_2)}\right]+mgl(2\cos{\theta_1}+\theta_2)
\)
と書けます。あとはラグランジュの運動方程式を当てはめて計算します。

3, 新たな座標系で運動方程式を導く
保存力場中でのラグランジュの運動方程式は
\(
\begin{eqnarray}
\left\{
\begin{aligned}
\frac{d}{dt}\left(\frac{\partial L}{\partial\dot{\theta}_1}\right)-\frac{\partial L}{\partial\theta_1}&=0 \\
\frac{d}{dt}\left(\frac{\partial L}{\partial\dot{\theta}_2}\right)-\frac{\partial L}{\partial\theta_2}&=0
\end{aligned}
\right.
\end{eqnarray}
\)
なので、代入し、\(\ddot{\theta}_1,\ddot{\theta}_2\)に関する運動方程式にすれば
\(
\begin{eqnarray}
\left\{
\begin{aligned}
\ddot{\theta}_1&=\frac{1}{2}\left\{-\ddot{\theta}_2\cos{(\theta_1-\theta_2)}-\dot{\theta}_2^2\sin{(\theta_1-\theta_2)}-2\frac{g}{l}\sin{\theta_1} \right\} \\
\ddot{\theta}_2&=-\ddot{\theta}_1\cos{(\theta_1-\theta_2)}+\dot{\theta}_1^2\sin{(\theta_1-\theta_2)}-\frac{g}{l}\sin{\theta_2}
\end{aligned}
\right.
\end{eqnarray}
\)
となります。これは非線形の2階連立微分方程式です。カオスです。解けません。
数値計算で解きます。用いるのは4次ルンゲ・クッタ法です。

実際に解いてみると、こうなります。

2重振り子の実験もされています。

[adsense2]

カオスが現れる場合、初期値鋭敏性という性質があります。
初期値がほんのちょっと変わるだけでその後の時間発展の様子が大きく変わる性質です。
ちょっと値を変えると上の動画とはまるで違う動きになります。
この動画は上の動画よりも中心に近いほうの粒子の初期速度が4%違います

機械精度の違いでさえも十分な時間発展の後には大きな違いが現れます。複雑であると同時に面白い現象です。

gnuplotで画像出力

gnuplotで画像を出力したいことがあります。
その方法を書きたいと思います。

想定する環境

  • gnuplot: Version 4.6 patchlevel 4
  • Imagemagick(convertコマンドを使用)

がインストールされている状況を考えます。特に、gnuplotではdofor構文を用いるため少なくともVersion4.6以上出ないと本稿に沿った画像出力はできませんのでご了承ください。

画像出力について

画像出力するための方法として、epsファイルを出力してからjpgやpng形式に変換します(僕がそれしか知らないためです)。

基本的には
ターミナル(x11とかwxtとかaquaとかpostscriptとかgifとか…)の指定→出力ファイル名の指定→グラフをファイルに出力
という流れになります。
手入力で画像出力する場合はgnuplotで

set terminal postscript eps enhanced color
set output "aaa.eps"
replot

とすれば今出力されているグラフの画像が”aaa.eps”というファイル名で得られます。
もしもjpgやpng形式に変換したければImagemagickのconvertコマンドを用いて、

convert -density 300x300 aaa.eps aaa.jpg

とすればjpg形式が得られます。jpgをpngに書き換えればpngが得られます。

しかし、何枚も画像がほしい場合には不向きです。スクリプトを書きましょう。
“out.plt”というファイルを作り、その中にこう書きます。

fname="$0"
set terminal postscript eps solid enhanced color
set output fname.".eps"
replot
set out
set terminal wxt enhanced

# 1. PNG with background
cv=sprintf("convert -density 300x300 %s.eps %s.png",fname,fname)

# 2. PNG with background
#cv=sprintf("convert -density 150x150 %s.eps %s.png",fname,fname)

# 3. PNG with background white
#cv=sprintf("convert -density 150x150 -background white -flatten -alpha off %s.eps %s.png",fname,fname)

system cv

このファイルをgnuplotのcallコマンドを用いて、

call "out.plt" "bbb"

として呼び出します。callは引数を指定できるコマンドです。
すると”bbb.eps”と”bbb.png”ファイルが作成されます。便利だね!

縮約

縮約はアインシュタインが一般相対性理論を説明するために最初に導入した数式のお約束です。
これが考え出されたのはアインシュタインが

和記号を書くのに飽きた

からであると予想されます。
このお約束なんですが,大学の物理科ではあまり触れられない癖に教授たちは常識のように語るので,
身につけておいて損はないですね。

さて,あるテンソル量\(A_{i}\)と\(B_{i}\)があったとしましょう。
この二つの縮約をとるということは,同じテンソル同士の和を取る,という意味になります。
アインシュタインの記法に則れば次のように書きます。

\(
A_{i}B^{i}=\sum_{i}{A_{i}B_{i}}
\)

ようするに同じ添え字が右下と右上に来たら和をとりましょうという約束です。
上の式は内積を表していることがわかります。
これを使うと外積の\(i\)成分も次のようになります。

\(
(A \times B)_i= \epsilon_{i,j,k}A^{j}B^{k}
\)
\(
\epsilon_{i,j,k}= 1 ~~ {\rm at} ~~ i,j,k=(1,2,3),(2,3,1),(3,1,2)
\)
\(
\epsilon_{i,j,k}= -1 ~~ {\rm at} ~~ i,j,k=(1,3,2),(3,2,1),(2,1,3)
\)
\(
\epsilon_{i,j,k}= 0 ~~ {\rm at} ~~ i,j,k=(otherwise)
\)

この\(\epsilon_{i,j,k}\)をレヴィチビタの完全反対称テンソルといいます。

モンテカルロ法

モンテカルロ法は,乱数を使って系の平衡状態にせまる一つの方法です。

ここでは二次元正方格子のイジング模型について紹介します。
イジング模型は隣り合う格子のスピン同士の相関エネルギーのみを考慮した模型です。
ハミルトニアンは以下のようになります。

\(
\displaystyle
H=-J \Sigma_{< i , j >}{s_{i}s_{j}}
\)
 
あるスピン状態を持つ系の状態Aについて,そのハミルトニアンを\(H_A\)とすると,
ボルツマン因子と分配関数\(Z\)を用いて状態の実現確率を計算することが出来ます。

\(
P(A)=\frac{\exp{(-\frac{H_A}{k_B T})}}{Z}
\)
 
Aからランダムに一つの格子におけるスピンを一つだけフリップした状態を状態Bとして,
状態Aと状態Bの実現確率比を乱数と比較します。
\(P(B)/P(A)>{\rm random}\)ならば状態Bを採用し,そうでなければ状態Aを採用します。
ここまでの流れを1mcs(モンテカルロステップ)といい,これを何度も繰り返すことで平衡状態に近づきます。
以下の二つのgifは,相転移温度より高い温度と低い温度で,
モンテカルロステップにより二次元正方格子のスピンがどのように振る舞うかをみているアニメーションです。
黒がダウンスピン,白がアップスピンです。


温度T=10eVのスピンの挙動
温度T=10eVのスピンの挙動

こちらは温度T=10eVで,自発磁化を持たない事が分かります。


温度T=0.001eVのスピンの挙動
温度T=0.001eVのスピンの挙動

こちらは温度T=0.001eVで,スピンが揃って自発磁化を持つ事が分かります。

画像をクリックするとmcsスタート!!

スピンのふるまいを追うと,自発磁化を持つ温度と持たない温度でスピンのふるまいが違う事が分かります。
二次元のイジング模型はオンサーガーにより厳密に解かれている(相転移温度T=2.26eV at J=1)ので,
それを踏まえて下のプログラムで遊んでみるといいですよ。

これのプログラムソースコード(保証は一切しません!!)
*注意 このプログラムを実行する前にising_anというディレクトリを作ってください。
*追記 2015 2/13 磁化とエネルギーを計算出来るようにしました。
計算方法は重み付き選択法というのを使ってます。

\(
\langle A \rangle = \frac{1}{T} \sum_{t=t_0+1}^{t_0 + T}A_{\alpha_{t}}
\)
 
十分平衡状態になっている\(t_0\)以降のモンテカルロステップで,物理量の統計平均を取る方法です。
平衡以前の物理量を計算しないmcsをからまわしと言います。

格子点51×51で,\(J=2\)で計算しました。
スピンの大きさを1として,相互作用をダブルカウントしているので注意しましょう。
即ち\(J=1\)の結果が欲しければプログラム上ではintJ=0.125にすればよいかも。

**********parameter**********
i_max::x方向の格子点の数
j_max::y方向の格子点の数
mcs_max::最大モンテカルロステップ数
mcs2step::mcs2step以降のmcsを平衡状態と見なし,物理量を計算します
intj::\(J/4\)
temp::温度
mov_max::mcs_max/movmax毎のスピン状態をgifアニメーションで見れるように区切ります。
**********parameter**********

単位格子あたりのエネルギーの温度依存性グラフを添付します。

ene

相転移温度付近でエネルギーの変化が大きい事が分かります。

単位格子あたりの磁化の絶対値の温度依存性グラフを添付します。

mag

相転移温度付近で自発磁化が発生し始めている事が分かります。

gnuplotを用いて,格子点上のスピンがどのように平衡状態へ近づくかをダイナミックに確認する事が出来ます。

gnuplot anime.gnuplot
 
を実行するとアニメが始まります。

gnuplot animegif.gnuplot
 
でgifアニメが作成出来ます。

!————————————————
! title = ising2d.f90
! developer = Amesyabody
! released = 2015 2/13
! programming language = fortran
!
!explanation
!2 dimentional square lattice Ising model's dynamical magnetism
!calculated by Monte Carlo method.
!————————————————
   
program main
    implicit none
    integer(4)::i,j,i_max,j_max,mcs,mcs_max,mov,mov_max,mcs2step
    parameter(i_max=50,j_max=50,mcs_max=1000000,&
    mcs2step=500000,mov_max=100)
    real(8)::intj,temp
    parameter(intj=0.25d0,temp=0.1d0)
    real(8)::spin(0:i_max,0:j_max)
    real(8)::ratio,ham1,ham2,hamres,dum,ene
    integer(4)::px,mx,py,my,rndx,rndy
    character(10)::movc,movc2

    real(8)::mag,mag2,magres
   
    !-----make random-----
    integer(4),allocatable::seed(:)
    integer cnt,seedsize
    real(8)::rnd,rndi,rndj,rnds
   
    call random_seed(size=seedsize)
    allocate(seed(seedsize))
    write(*,*)seedsize

    dum=rndmf(0)

    call system_clock(count=cnt)
    seed=cnt!+idnint(dum*1000000)
!            dum=rndmf(seed(0))
!            seed=idnint(seed*dum)
    call random_seed(put=seed)

    !------------------

mag2=0.0d0
ene=0.0d0

do mov=0,mov_max

write(movc,'(i3)')mov

open(100+mov,file="ising_an/spin"//trim(adjustl(movc))//".dat"&
,status="unknown")

end do

open(5000,file="ising_an/anime.gnuplot",status="unknown")
open(5001,file="ising_an/animegif.gnuplot",status="unknown")

    !-----primitive spin set-----
   
    do i=0,i_max
        do j=0,j_max
            call random_number(rnd)
!           write(*,*)rnd
            if(rnd>0.5d0)then
                spin(i,j)=1.0d0
            else
                spin(i,j)=-1.0d0
            end if
!           write(*,*)spin(i,j)!,cnt,seed(1)

        write(100,'(i3,i3,f5.1)')i,j,spin(i,j)

        end do

        write(100,*)" "

    end do
   
    !----------------------------
   
    do mcs=1,mcs_max
   
        do i=0,i_max
       
            do j=0,j_max
           
                px=i+1
                if(i==i_max)px=px-i_max-1
                py=j+1
                if(j==j_max)py=py-j_max-1
                mx=i-1
                if(i==0)mx=mx+i_max+1
                my=j-1
                if(j==0)my=my+j_max+1

                ham1=ham1-intj*spin(i,j)*spin(px,j)&
                -intj*spin(i,j)*spin(i,py)&
                -intj*spin(mx,j)*spin(i,j)&
                -intj*spin(i,my)*spin(i,j)

                do mov=1,mov_max-1

                    if(mcs==mov*mcs_max/mov_max)&
                    write(100+mov,'(i3,i3,f5.1)')i,j,spin(i,j)

                end do

                if(mcs>mcs2step)then

                    mag2=mag2+spin(i,j)

                end if
               
            end do

            do mov=1,mov_max-1

                if(mcs==mov*mcs_max/mov_max)write(100+mov,*)" "

            end do

        end do

        if(mcs>mcs2step)then

            magres=magres+abs(mag2)
            mag2=0.0d0
            ene=ene+ham1

        end if

        call random_number(rndi)
        call random_number(rndj)
        rndx=idnint(rndi*i_max)
        rndy=idnint(rndj*j_max)

!       write(*,*)rnd
!write(*,*)spin(rndx,rndy)
        spin(rndx,rndy)=-spin(rndx,rndy)
!write(*,*)spin(rndx,rndy)
        do i=0,i_max
       
            do j=0,j_max
           
                px=i+1
                if(i==i_max)px=px-i_max-1
                py=j+1
                if(j==j_max)py=py-j_max-1

                mx=i-1
                if(i==0)mx=mx+i_max+1
                my=j-1
                if(j==0)my=my+j_max+1

                ham2=ham2-intj*spin(i,j)*spin(px,j)&
                -intj*spin(i,j)*spin(i,py)&
                -intj*spin(mx,j)*spin(i,j)&
                -intj*spin(i,my)*spin(i,j)
   
            end do
           
        end do
       
        ratio=exp((-ham2+ham1)/temp)
!       write(*,*)ham1,ham2,ratio
       
        ham1=0.0d0
        ham2=0.0d0
       
        call random_number(rnds)

        if(ham1 > ham2)then!metroporis
        goto 200
        end if
       
        if(rnds>ratio)then
        spin(rndx,rndy)=-spin(rndx,rndy)
        end if

200     continue

    end do
   
    do i=0,i_max
        do j=0,j_max
       
        mag=mag+spin(i,j)
       
        px=i+1
        if(i==i_max)px=px-i_max-1
        py=j+1
        if(j==j_max)py=py-j_max-1
        mx=i-1
        if(i==0)mx=mx+i_max+1
        my=j-1
        if(j==0)my=my+j_max+1

        hamres=hamres-intj*spin(i,j)*spin(px,j)&
        -intj*spin(i,j)*spin(i,py)&
        -intj*spin(mx,j)*spin(i,j)&
        -intj*spin(i,my)*spin(i,j)

        write(100+mov_max,'(i3,i3,f5.1)')i,j,spin(i,j)
       
        end do

        write(100+mov_max,*)" "

    end do
   
    mag=mag/(dble(i_max*j_max))
   
!   write(*,'(3f10.5,f15.5)')temp,intj,mag,hamres

    ene=ene/(dble((mcs_max-mcs2step)*(i_max+1)*(j_max+1)))
    magres=magres/(dble((mcs_max-mcs2step)*(i_max+1)*(j_max+1)))

    write(*,'(A,f10.5,A,f10.5)')"energy=",ene,"      magnet=",magres

    write(5000,*)"set size sq"
    write(5000,*)"set pm3d map"
    write(5000,*)"set palette rgbformulae 21,22,23"
    write(5000,*)"set xlabel 'x'"
    write(5000,*)"set ylabel 'y'"
    write(5000,*)"set zra[-1:1]"
    write(5000,*)"set zlabel 'spin'"
    write(5000,*)"splot 'spin0.dat'"
    write(5000,*)"pause 0.1"
    do mov=1,mov_max
    write(movc,'(i3)')mov
    write(5000,*)"splot 'spin"//trim(adjustl(movc))//".dat'"
    write(5000,*)"pause 0.1"

    end do

    write(5001,*)"set tics font 'Times,15'"
    write(5001,*)"set xlabel font 'Times,15'"
    write(5001,*)"set ylabel font 'Times,15'"
    write(5001,*)"set size sq"
    write(5001,*)"set pm3d map"
    write(5001,*)"set palette rgbformulae 21,22,23"
    write(5001,*)"set xlabel 'x'"
    write(5001,*)"set ylabel 'y'"
    write(5001,*)"set zra[-1:1]"
    write(5001,*)"set cbrange [-1:1]"
    write(5001,*)"set zlabel 'spin'"
    write(5001,*)"set term gif animate optimize"
    write(5001,*)"set output 'is_t10.gif'"
    write(5001,*)"do for[i=0:100] {"
    write(5001,*)'file = sprintf("spin%01d.dat", i)'
    write(movc2,'(i10)')mcs_max/mov_max
    write(5001,*)'time = sprintf("s=%d[mcs]",i*'//trim(adjustl(movc2))//')'
    write(5001,*)"set title time"
    write(5001,*)"splot file"
    write(5001,*)"}"

    contains

real(8) function rndmf(seeds)
implicit none

integer(4)::a,m,q,p,n,ndiv,j,k,seeds
real(8)::rm,rmax
parameter(a=16807,m=2147483647,rm=1.0/m)
parameter(q=127773,p=2836,n=32,ndiv=1+(m-1)/n)
parameter(rmax=1.0-1.2e-7)
integer(4)r(n),r0,r1

if (seeds .ne. 0)then
r1=abs(seeds)
do j=n+8,1,-1
k=r1/q
r1=a*(r1-k*q)-p*k
if(r1 .lt. 0)r1=r1+m
if(j .le. n)r(j)=r1
end do
r0=r(1)
end if


k=r1/q
r1=a*(r1-k*q) -p*k
if(r1 .lt. 0)r1=r1+m
j=1+r0/ndiv
r0=r(j)
r(j)=r1
rndmf=min(rm*r0,rmax)

end function
   
end program