R版KETpic

2011年11月13日 (日)

R版KETpicタートルグラフィックスによる図形作成(作成中)

R版KETpicに組み込んだタートルグラフィックスによる描画

本記事の一部は「第34回CASTeXセミナー」で紹介する予定です.
日時:平成23年12月25日(日)15:00~15:55

※使用するketpicパッケージなどについては,本ページの
『対数正規確率紙および正規確率紙の作成とプロットデータ例』の項を見て下さい.
ketpicパッケージ等は更新することがありますので注意して下さい

■作図例  左から1),2),3)

  Hoshi2      Pentagon2     Re6polygon


図中の三角(45°三角定規のかたち)はタートル(カメ)で,作図完了時のカメの位置と方向を示している.ここでは,黒色で描いたものを表示.

■正n多角形の一つの内角の大きさ(60分法,度数法,45°)は次式で表される.

Sei_n_kakkei_naik_ookisa

この式を用いて,正n多角形を描くコードを作成する.
・タートルの初期の向きは東方向(画面に向かって右方向)を向いている.
・使用コマンドは授業で説明.

■R Graphics画面(座標範囲は任意に設定可)

Tg_coord2   St_gamen1

■タートルグラフィックスによる作図により,条件判定文(if文),繰り返し文(for文,while文)を学習する.

■コード例
●作図例 1
rm(list=ls(all=TRUE))#ワークスペース上の全てのオブジェクトを消去
setwd("C:/rwork")
load("ketpic110824.RData")#準備中
load("tgpack111105.RData")#準備中
Tg_init()
WindispT()
NewKame(0,0)
SetHeading(72)
s<- 3
kaku <- 144
for (i  in 1:5) {
  Forwd(s)
  TurnLeft(-kaku)
}
TgDraw(color="red")#Rの描画ウインドウで赤色で表示

●作図例 2),3)
rm(list=ls(all=TRUE))#ワークスペース上の全てのオブジェクトを消去
setwd("C:/rwork")
load("ketpic110824.RData")#準備中
load("tgpack111105.RData")#準備中
Tg_init()
WindispT()
NewKame(0,0)
hen<- 2
n<- 5     # 作図例(2) : n<- 5 ,  作図例(3) : n<- 6
naikaku<- 180*(n-2)/n
for (i in 1:n){         #for文
  Forwd(hen);TurnLeft(180-naikaku)
}
TgDraw(color="black")

■うずまき
Uzu_shikaku
 

rm(list=ls(all=TRUE))#ワークスペース上の全てのオブジェクトを消去
setwd("C:/rwork")
load("ketpic110824.RData")#準備中
load("tgpack111106.RData")#準備中
Tg_init()
WindispT()
NewKame(0,0)
hen<- 0.5
step<- 0.2
for (i in 0:40){
  TurnLeft(90);Forwd(hen);
  hen<- hen+step
}
TgDraw()
############以下,図ファイルの作成(Texにinputする図ファイル)
Openfile("zu_uzu_shikaku.tex")
Beginpicture("1cm")
  TgDrwline()
Endpicture(0)
Closefile()
#############
■再帰呼び出しによるコード例
rm(list=ls(all=TRUE))#ワークスペース上の全てのオブジェクトを消去す
setwd("C:/rwork")
load("ketpic110824.RData")#準備中
load("tgpack111106.RData")#準備中

uzu<- function(hen,step){
  if (hen>10) return()
  else {
    TurnLeft(90);Forwd(hen)
    return(Recall(hen+step,step))  #return(uzu(hen+step,step))
  }
}

Tg_init()
WindispT()
NewKame(0,0)
hen<- 0.5
step<- 0.2
uzu(hen,step)
TgDraw(color="red",width=2)

■花火

Hanabi

rm(list=ls(all=TRUE))#ワークスペース上の全てのオブジェクトを消去す
setwd("C:/rwork")
load("ketpic110824.RData")#準備中
load("tgpack111107.RData")#準備中
#hanabi
Tg_init()
WindispT()
NewKame(0,0)
s<- 5
kaku <- 36
for (i in 1:10) {
  TurnLeft(kaku)
  Forwd(s);TgDraw(color="red",width=2);
  Puttext(Where(),i,pos=3)
  Back(s)
  Sleep(1) #1秒休む(ゆっくり描画)
}

■斜め投射運動(タートルグラフィックスによる)
初速度(v0)=4,7,8,9m/s ,投げる角度theta=45°の軌跡の様子
0.1秒間隔で1.4秒間の軌跡変化を示す.小さい三角はタートルを表示.

Nanamenage

############
●コード例(初速度(v0)=4,7m/s 投げる角度theta=45°の場合)
rm(list=ls(all=TRUE))#ワークスペース上の全てのオブジェクトを消去
setwd("C:/rwork")
load("ketpic110824.RData")
load("tgpack111108.RData")
Tg_init(c(0,10),c(-3,3))
WindispT()
Kameinit(s)#カメのサイズを指定
NewKame(0,0)
s<- 0.4
time<- 14
G<- 9.8 #m/s^2
theta0<- 45 #投げる角度
v0<- 4 #初速度m/s
theta<- theta0*pi/180
for (i in 1:time){
  t<- i/10
  x<- v0*cos(theta)*t;
  y<- v0*sin(theta)*t-1/2*G*t*t
  LineTo(x,y)
  TgDraw(color="blue",width=2)
}
#########
Tg_init(c(0,10),c(-3,3))
Kameinit(s)
NewKame(0,0)
v0<- 7 #初速度m/sm/s
theta0<- 45*pi/180
for (i in 1:time){
  t<- i/10
  x<- v0*cos(theta)*t;
  y<- v0*sin(theta)*t-1/2*G*t*t
  LineTo(x,y)
  TgDraw(color="blue",width=2)
}
Puttext(c(5,2.8),"斜め投射運動")
Puttext(c(5,2.4),paste("初速度(v0)=",v0,"m/s"," 投げる角度theta=",round(theta0,digits=1),"°",sep=""))

■時計 

Tgclock

■時計のコード例 1    1分間隔で表示.1時間まで表示すると思います.
rm(list=ls(all=TRUE))#ワークスペース上の全てのオブジェクトを消去
setwd("C:/rwork")
load("ketpic110824.RData")
load("tgpack111108.RData")
kaku30<- 30
runtime_min<- 1
stime<- as.POSIXlt(Sys.time())
runtime_hours<- 0
while (runtime_hours<=1){#1hour
  time<- as.POSIXlt(Sys.time())
  hour<- time$hour %% 12
  min<- time$min
  if(runtime_min>=1){
    stime_min<- as.POSIXlt(Sys.time())
    Lkaku<- min*360/60 #long_hand(minute_hand)
    Skaku<- hour*kaku30 +  1/12*Lkaku)#short_hand(hour_hand)
   
    Tg_init(c(-5,5),c(-5,5))
    NewKame(0,0)
    WindispT()
    SetHeading(90)
    TurnRight(Skaku)##hour_hand
    Forwd(1);TgDraw(color="blue",width=2);Back(1)
 
    Tg_init(c(-5,5),c(-5,5))
    NewKame(0,0)
    SetHeading(90)
    TurnRight(Lkaku)#minute_hand
    Forwd(2);TgDraw(color="red",width=2);Back(2)
  }
  runtime_hours<- difftime(time,stime,unit="hours")
  runtime_min<- difftime(time,stime_min,unit="min")

  en<- Circledata(c(0,0),2.5)
  WindispT(en,new=TRUE)
  memori<- c(1,2,3,4,5,6,7,8,9,10,11,12)
  r<- 2.2
  theta<- 60*pi/180
  for (i in 1:12){
    Puttext(c(r*cos(theta),r*sin(theta)),memori[i])
    theta<- theta-30*pi/180
  }
}

■時計のコード例2   秒間隔で表示.2時間まで表示すると思います.
rm(list=ls(all=TRUE))#ワークスペース上の全てのオブジェクトを消去
setwd("C:/rwork")
load("ketpic110824.RData")
load("tgpack111108.RData")

kaku30<- 30
runtime_min<- 1
stime<- as.POSIXlt(Sys.time())
runtime_hours<- 0
while (runtime_hours<=2){#2時間表示
  time<- as.POSIXlt(Sys.time())
  hour<- time$hour %% 12
  min<- time$min
  sec<- time$sec

  Lkaku<- min*6 + sec*6/60 #long_hand(minute_hand)
  Skaku<- hour*kaku30 + 1/12*Lkaku#*short_hand(hour_hand)
  Tg_init(c(-5,5),c(-5,5))
  NewKame(0,0)
  WindispT()
  SetHeading(90)
  TurnRight(Skaku)#hour_hand
  Forwd(1);TgDraw(color="blue",width=2);Back(1)
 
  Tg_init(c(-5,5),c(-5,5))
  NewKame(0,0)
  SetHeading(90)
  TurnRight(Lkaku)#minute_hand
  Forwd(2);TgDraw(color="red",width=2);Back(2)

  en<- Circledata(c(0,0),2.5)
  WindispT(en,new=TRUE)
  memori<- c(1,2,3,4,5,6,7,8,9,10,11,12)
  r<- 2.2
  theta<- 60*pi/180
  for (i in 1:12){
    Puttext(c(r*cos(theta),r*sin(theta)),memori[i])
    theta<- theta-30*pi/180
  }

  time<- as.POSIXlt(Sys.time())
  runtime_hours<- difftime(time,stime,unit="hours")#表示時間の取得
}

■カラーの箱

Color_box

rm(list=ls(all=TRUE))#ワークスペース上の全てのオブジェクトを消去
setwd("C:/rwork")
load("ketpic110824.RData")
load("tgpack111112.RData")

shikaku2<- function(...,dx,dy){#四角を描く関数.四角の中心の座標を与える(x,y)
  L<- list(...)
  if (length(L)==1){
    Tmp<- L[[1]]
    x<- Tmp[1];y<- Tmp[2]
  }
  else{
    x<- L[[1]];y<- L[[2]]
  }
  p<- c(x-dx,y-dy)
  NewKame(p)
  Forwd(2*dx);TurnLeft(90)
  Forwd(2*dy);TurnLeft(90)
  Forwd(2*dx);TurnLeft(90)
  Forwd(2*dy)
}

Tg_init()
hen<- 1
shikaku2(c(2,2),dx=hen,dy=hen)
shikaku2(c(2,-2),dx=hen,dy=hen)
shikaku2(c(-2,2),dx=hen,dy=hen)
shikaku2(c(-2,-2),dx=hen,dy=hen)

KameMat()
for(i in 1:4){
  col<- rainbow(12)[i]
  TgPaint(TgPaintdata(i),col=col,border=col)
}

■カラーパレット

Color_palete  Heat_colors   Cm_colors  Topo_colors  Terrain_colors

rm(list=ls(all=TRUE))#ワークスペース上の全てのオブジェクトを消去
setwd("C:/rwork")
load("ketpic110824.RData")
load("tgpack111112.RData")

shikaku2<- function(...,dx,dy){#四角を描く関数.四角の中心の座標を与える(x,y)
  L<- list(...)
  if (length(L)==1){
    Tmp<- L[[1]]
    x<- Tmp[1];y<- Tmp[2]
  }
  else{
    x<- L[[1]];y<- L[[2]]
  }
  p<- c(x-dx,y-dy)
  NewKame(p)
  Forwd(2*dx);TurnLeft(90)
  Forwd(2*dy);TurnLeft(90)
  Forwd(2*dx);TurnLeft(90)
  Forwd(2*dy)
}

Tg_init()
r<- 4
hen<- 0.5
theta<- 0
for(i in 1:20){
  P<- r*c(cos(theta),sin(theta))
  shikaku2(P,dx=hen,dy=hen)
  theta<- theta+18*pi/180
}

KameMat()
for(i in 1:20){
  col<- rainbow(20)[i] #heat.colors(20)[i] cm.colors(20)[i] topo.colors(20)[i] terrain.colors(20)[i]
  TgPaint(TgPaintdata(i),col=col,border=col)
}

■RGB Color,HSV Color

Color_rgb   Color_hsv   

Setwindow(c(-1,1),c(-1,1))
WindispT()
for (i in 1:1000){
  Pt<- runif(2,min=-1,max=1)
  r<- runif(1,min=0,max=1) #h
  g<- runif(1,min=0,max=1) #s
  b<- runif(1,min=0,max=1) #v
  colrgb<- rgb(r,g,b)          #colhsv<- hsv(h,s,v)
  WindispT(Pt,color=colrgb,pch=20,cex=3,type="p",new=TRUE)
}

■円周率πの近似計算例
「The R Tips 」舟尾暢男著 p84  2005.3.1 九天社 を利用させていただきました.

Montecarlo     Montecarlo2

setwd("C:/rwork")
load("ketpic110824.RData")
load("tgpack111117.RData")

pi.montecarlo <- function(n) {
   if(n>50000){
      print("回数は50,000以下にして下さい")
      return()
   }
   Writeln("回数は",n,"回です")
   count <- 0         
   for (i in 1:n) {
     suna <- runif(2)   
     G<- Pointdata(suna)
     WindispT(G,color=gray(10/12),pch=20,cex=0.05,type="p",new=TRUE)
     if (sqrt(suna[1]^2 + suna[2]^2) < 1) {
       count <- count + 1
     }
   }
   Writeln("円周率π=",4 * count / n)
   return(4 * count / n)
}

Setwindow(c(-1,1),c(-1,1))
WindispT()
n<- 50000
pai<- pi.montecarlo(n)
en<- Circledata(c(0,0),1, "R=c(0,pi/2)")
WindispT(en,color="red",new=TRUE)
Puttext(c(-0.5,0.6),"関数pi.montecarloの")
Puttext(c(-0.5,0.5),paste("実行回数は",n,"回です",sep=""))
Puttext(c(-0.5,0.25),paste("円周率π=",pai,sep=""))

■出力文について(学習用の自作関数として限定的な活用として使用します)
Write文,Writeln文:改行なし,改行あり
Writef文,Writelnf:改行なし,改行あり.書式設定付
【利用例】
・斜め投射運動の軌跡の座標値算出
###########
g<- 9.8;
v0<- 9;
theta0<- 45;
theta<- 45*pi/180;
Writelnf("%3s %5s %5s","秒","x","y")
for (i in 0:15){
  t<-  i/10
  x<-  v0*cos(theta)*t
  y<-  v0*sin(theta)*t-1/2*g*t*t
  t<- sprintf("%3.1f",t)
  x<- sprintf("%5.2f",x);
  y<- sprintf("%5.2f",y);
  Writeln(t, "," , x , "," , y)
}
###########
参考 上記コードの内,
  t<- sprintf("%3.1f",t)
  x<- sprintf("%5.2f",x);
  y<- sprintf("%5.2f",y);
  Writeln(t, "," , x , "," , y)
の4行は,下記のWritelnf文1行に置き換えることができる.
  Writelnf("%3.1f%s%5.2f%s%5.2f",t,",",x,",",y)
###########
g<- 9.8;
v0<- 9;
theta0<- 45;
theta<- 45*pi/180;
Writelnf("%3s %5s %5s","秒","x","y")
for (i in 0:15){
  t<-  i/10
  x<-  v0*cos(theta)*t
  y<-  v0*sin(theta)*t-1/2*g*t*t
  Writelnf("%3.1f%s%5.2f%s%5.2f",t,",",x,",",y)
}
###########
【実行結果】
Houbutu_kisekichi

#####
■対数正規確率紙および正規確率紙の作成とプロットデータ例
【下上段右図】『環境工学系のための数学』 滝沢智 著 数理工学社20041210
p27 例題1.14の作図例

Lognormaldist  Normaldist

Lognormaldist_b_2
各軸の格子線の範囲と目盛りの値ラベルは任意に指定.
データの最小値と最大値が表示されるので,この値を参考に指定.
min=  0.972421758592559  Cumulative Percent= 1
max=  200.950838178908  Cumulative Percent= 99

縦軸の確率線:c( 0.1, 1, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 99, 99.9)の15本

■ソースコード例(わかりにくい箇所や不備な箇所があると思いますが参考に掲載しました)
ketpic110824.RData
tgpack111205.RData
以下からDLできます.(ユーザ名:  knctknct   パスワード: ecoeco )
http://kiyomidai.cocolog-nifty.com/blog/2011/12/post-3619.html

本ページのプログラムを『ketpic110824.RData』と『tgpack111205.RData』
を利用するものに修正して下記に掲載しました.
(ユーザ名:  knctknct   パスワード: ecoeco )

http://kiyomidai.cocolog-nifty.com/blog

●対数正規確率紙
作成された図ファイル(zu_lognormal.tex)はtex文書に挿入して利用.
texについてはHPや図書館で学習してください.
\input{zu_lognormal.tex}
 「lnpp111207_b.r」をダウンロード

●正規確率紙
図ファイル(zu_normal.tex)
「npp111207.r」をダウンロード

下記のページを参考にさせていただきました.
・対数正規確率紙
http://aoki2.si.gunma-u.ac.jp/R/lnpp.html
・正規確率紙 
http://aoki2.si.gunma-u.ac.jp/R/npp2.html
・確率紙プロット
http://www.okada.jp.org/RWiki/?%B3%CE%CE%A8%BB%E6%A5%D7%A5%ED%A5%C3%A5%C8

2011年10月26日 (水)

R版KETpicでフローチャートを描く/作表例(作成中)

■R版KETpicでフローチャート(Pascal資料用)を描いてみよう.
 プログラムによる作図のため,修正と再利用が容易になると思います.

■作図例A  
  右図:繰り返し文(定回反復,for文)のフローチャート
  左の逐次構想と右の反復構造は同じ内容である.

     Flowchart1          Flowchart2

■作図例B(111029追加)

 繰り返し文(不定回反復,while文とrepeat文)のフローチャート
         【while文】                     【repeat文】

  While         Repeat_until

・ひし形部に記述する条件は何を判定するのかを明確にする.
単に条件としないで,継続条件または終了条件とする.

■描き方の方針
1) 四角形,ひし形を配置する.
   ・四角形のデータには,図心と4辺の中点の座標値をもたせる.
   ・ひし形のデータには,図心と4つの頂点の座標値をもたせる.
2) 線は四角形やひし形の図形上を通して引く.
3) 四角形やひし形の図形部を塗りつぶし(色設定も可)不必要な線を消す.
4) 必要な座標位置に矢印を配置する.
5) 文字列を適切な座標位置に配置する.

■パッケージのDL(準備中)
ketpic110824.RData
tgpack111029.RData

■作図例Aの左図のコード例

#Tg_packの利用 図ファイルzu_flowchart1.tex
#NewKame,SetPosition,SetHeading,LineTo,Tg_listなどはタートルグラフィックスのコマンド
#図は座標(0,0)を中心に描き始める

rm(list=ls(all=TRUE))#ワークスペース上の全てのオブジェクトを消去する
setwd("C:/rwork")
load("ketpic110824.RData")#更新作業中
load("tgpack111029.RData")#更新作業中
setwd("C:/rwork")

Setwindow(c(-1.5,3),c(-5,1))
Openfile("zu_flowchart1.tex")
Beginpicture("1cm")

#タートルグラフィックスの初期化
Tg_init()

#左の四角形を5個描く(s:=0など)
box0all<- list()
dy<- 0
for (i in 0:4){
  box0<- Boxdata(c(0,-dy),1,0.25)
  box0all<- c(box0all,box0)
  dy<- dy+1
  Drwline(box0[[1]])
}

#右の四角形を5個描く(数字が記入された箱)
box1all<- list()
dy<- 0
for (i in 0:4){
  box1<- Boxdata(c(2,-dy),0.5,0.25)
  box1all<- c(box1all,box1)
  dy<- dy+1
  Drwline(box1[[1]])
}

#亀の生成
NewKame(0,0)
SetPosition(0,1);SetHeading(-90)

LineTo( box0all[[30]][1],box0all[[30]][2]-1)#現在の位置から指定座標まで線を書く.
Tg_list<- Flattenlist(Tg_list)
Drwline(Tg_list)

#shodeで四角形の中の線を消す
for (i in 1:5) {
  BoxShdata(box0all[[6*i]],1,0.25,"white",1)
}

#矢印を配置
Arrowhead(Tg_list[[1]][2,],Tg_list[[1]])
for (i in 1:5){
  Arrowhead(matrix(box0all[[3+(i-1)*6]],ncol=2),Tg_list[[1]])
}

#大きな囲枠の四角形を描く
Waku<- Framedata(box0all[[18]],2,2)
Waku<- Translatedata(Waku,0.7,-0.5)#Wakuの平行移動(外枠の位置調整のため)
Drwline(Waku)

#左の四角形の中に文字列を配置 Fcletter文字列配置関数
dy<- 0
for (i in 0:4){
  Fcletter(c(0,dy),"c",paste("s:=",as.character(i),";",sep=""),"L")
  dy<- dy-1
}

#右の四角形の中に文字列を配置
dy<- 0
for (i in 0:4){
  Fcletter(c(2,dy),"c",as.character((i*(i+1)) %/% 2),"L")
  dy<- dy-1
}

#右の四角形の上部にsを配置
Fontsize("L")
Letter(box1all[[3]],"n2","s")
Fontsize("n")

Endpicture(0)
Closefile()

■作図例Aの右図のコード例

#Tg_packの利用 図ファイルzu_flowchart1.tex
#NewKame,SetPosition,SetHeading,LineTo,Tg_listなどはタートルグラフィックスのコマンド
#図は座標(0,0)を中心に描き始める

rm(list=ls(all=TRUE))#ワークスペース上の全てのオブジェクトを消去する
setwd("C:/rwork")
load("ketpic110824.RData")#更新作業中
load("tgpack111029.RData")#更新作業中
setwd("C:/rwork")

#Tg_packの利用  zu_flowchart2.tex
Setwindow(c(-2.5,3),c(-5,2))
Openfile("zu_flowchart2.tex")
Beginpicture("1cm")

#タートルグラフィックスの初期化
Tg_init()
#亀の生成
NewKame(0,0)

#ひし形と2つの四角図(box)の座標値を取得
dia0<- Diadata(c(0,0),1.5,0.7)
box10<- Boxdata(dia0[[5]]+c(0,-1),1,0.25)#ひし形の下のbox
box11<- Boxdata(dia0[[3]]+c(0,1),2.2,0.25)#ひし形の上のbox

#繰り返し部の線データを取得.タートルグラフィックスを利用
SetPosition(box11[[5]][1],box11[[5]][2]);SetHeading(-90)#亀のスタート準備
Forwd(4);TurnRight(90)#線1
Forwd(2);TurnRight(90)#線2
cyuten<- (box11[[5]]+dia0[[3]])/2
LineTo(WhereX(), cyuten[2]);TurnRight(90)#線3
LineTo(cyuten[1],cyuten[2])#線4

#ひし形から出るfalse部の線データを取得
SetPosition(dia0[[2]][1],dia0[[2]][2])
Forwd(1);TurnRight(90)#線5
Forwd(3)#線6

#線データを描く
Tg_list<- Flattenlist(Tg_list)
Drwline(Tg_list)

#ひし形とその下のboxの中の線を消す
BoxShdata(dia0[[5]]+c(0,-1),1,0.25,"cyan",0.5)#box中の線を消ため色cyanで塗る
DiaShdata(c(0,0),1.5,0.7,"white",1)#ひし形中の線を消すためにwhiteで塗る

#2つのboxとひし形を描く
Drwline(dia0[[1]],box10[[1]],box11[[1]])

#上のboxに文字列を配置
Fcletter(dia0[[3]]+c(0,1),"c","s:=0; //初期値の設定;","l")

#ひし形上に文字列を配置.2行の文字列の配置
Fcletter2(c(0,0),"n","継続条件","l","black","s","for i:=1 to 4","n","black")

#下のboxに文字列を配置
Fcletter(dia0[[5]]+c(0,-1),"c","s:=s+i;","L")

#矢印を配置
Arrowhead(Tg_list[[4]][2,],Tg_list[[4]])
Arrowhead(Tg_list[[6]][2,],Tg_list[[6]])

#フロー線に配置する文字列の設定
Fcletterrot((Tg_list[[3]][1,]+Tg_list[[3]][2,])/2,c(0,1),0,4,"繰り返し","s","red")
Fcletter(c(dia0[[2]][1],dia0[[2]][2]),"ne","FALSE","ss","red")
Fcletter(c(dia0[[5]][1],dia0[[5]][2]),"sw","TRUE","ss","red")

Endpicture(0)
Closefile()

■R版KETpicによる九九の表の作図例

Kuku

■上図のコード例

rm(list=ls(all=TRUE))#ワークスペース上の全てのオブジェクトを消去す
setwd("C:/rwork")
load("ketpic110824.RData")#更新作業中
load("tgpack111029.RData")#更新作業中
setwd("C:/rwork")

Setwindow(c(-5,55),c(-5,55))
WindispT()

tw<- 5
yw<- 5

tate<- as.list(rep(tw,10))
yoko<- as.list(rep(yw,10))

hyou<- Tabledata(c(-1,-1,5,5,5,5),tate,yoko)
#hyou<- Tabledata(c(-1,-1,5,5,5,0),tate,yoko)

G1<- hyou[[1]]# 表のPDデータ
syasen<- Diagcelldata(hyou,1,1)#左上のセル(1,1)の対角線PD syasen[[1]]

#セルの情報list(中心,横幅/2,縦幅/2)を返す
hyou_yoko_center<- Findcell(hyou,c(2,11),1) # 表の中心(行方向)j
hyou_tate_center<- Findcell(hyou,1,c(2,11)) # 表の中心(列方向)i

Tmp<- Findcell(hyou,1,1)
Pt0<- Tmp[[1]]+c(-Tmp[[2]],Tmp[[3]])#Pt0は表の左上隅の点

Pos11<- Findcell(hyou,1,1)
Pt11<- Pos11[[1]]#Pt11は表のセル(1,1)の斜線の中点

V<- yw/2*c(cos(pi/4),sin(pi/4))

Gs<- Listplot(Pt0,Pt0+V)#表から上に突き出た40度傾いた線のPD

#表から上に突き出た40度傾いた線のPD 11本
Gsall<- list()
for (i in 1:11){
  P<- Pt0+c((i-1)*tw,0)
  Gtmp<- Listplot(P,P+V)
  Gsall<- c(Gsall,list(Gtmp))
}

P1<- Ptend(Gsall[[1]])
Q1<- Ptend(Gsall[[11]])
Gyoko<- Listplot(P1,Q1)#表から上に突き出た部分の横線

#表から右側の突き出る線PD 11本
Pt2<- Ptstart(Gsall[[11]])
Gsall2<- list()
for (j in 1:11){
  P<- Pt2+c(0,-(j-1)*yw)
  Gtmp2<- Listplot(P,P+V)
  Gsall2<- c(Gsall2,list(Gtmp2))
}

P2<- Ptend(Gsall2[[1]])
Q2<- Ptend(Gsall2[[11]])
Gtate<- Listplot(P2,Q2)#表から右に突き出た部分の縦線

Tmp1<- Ptstart(Gsall[[1]])
Tmp2<- Ptstart(Gsall[[2]])
Tmp3<- Ptend(Gsall[[2]])
Tmp4<- Ptend(Gsall[[1]])
Gtmp<- Listplot(Tmp1,Tmp2,Tmp3,Tmp4,Tmp1)
Gh<- Hatchdata("i",list(Gtmp),45,5)#突き出た上面の左セルのハッチPD

Tmp1<- Ptstart(Gsall[[2]])
Tmp2<- Ptstart(Gsall[[11]])
Tmp3<- Ptend(Gsall[[11]])
Tmp4<- Ptend(Gsall[[2]])
Ue_men<- Listplot(Tmp1,Tmp2,Tmp3,Tmp4,Tmp1)#上に突き出た面のハッチ用PD
Gh_ue_men<- Hatchdata("i",list(Ue_men),45,5)#上に突き出た面のハッチ

Tmp1<- Ptstart(Gsall2[[1]])
Tmp2<- Ptstart(Gsall2[[2]])
Tmp3<- Ptend(Gsall2[[2]])
Tmp4<- Ptend(Gsall2[[1]])
Yoko_men<- Listplot(Tmp1,Tmp2,Tmp3,Tmp4,Tmp1)#横に突き出た面のうち,上1個のハッチ用PD
Gh_yoko_men<- Hatchdata("i",list(Yoko_men),45,5)#横に突き出た面のうち,上1個のハッチ用

Tmp<- Findcell(hyou,c(2,11),c(1,2))
YokoH<- Framedata(Tmp[[1]],Tmp[[2]],Tmp[[3]])#表1行の1から9のセルのハッチ用PD
YokoHach<- Hatchdata("i",list(YokoH),45,5)

Tmp<- Findcell(hyou,c(1,2),c(2,11))
TateH<- Framedata(Tmp[[1]],Tmp[[2]],Tmp[[3]])#表1列の1から9のセルのハッチ用PD
TateHack<- Hatchdata("i",list(TateH),45,5)#表1列の1から9のセルのハッチ

WindispT(G1,syasen[[1]],color="red",new=TRUE)
WindispT(Pt0,hyou_yoko_center[[1]],hyou_tate_center[[1]],color="blue",new=TRUE)
WindispT(Pt11,color="red",new=TRUE)
WindispT(Gs,color="blue",new=TRUE)
WindispT(Gsall,color="blue",new=TRUE)
WindispT(Gyoko,color="blue",new=TRUE)
WindispT(Gsall2,color="blue",new=TRUE)
WindispT(Gtate,color="blue",new=TRUE)
WindispT(Gh,color="blue",new=TRUE)
WindispT(Gh_ue_men,color="blue",new=TRUE)
WindispT(Gh_yoko_men,color="blue",new=TRUE)
WindispT(YokoHach,TateHack,color="blue",new=TRUE)

####
Openfile("zu_kuku_hyou.tex")
Beginpicture("1mm")
Setcolor("red")

Shade(Ue_men,0.2)
Shade(Yoko_men,0.2)

Shade(YokoH,0.2)
Shade(TateH,0.2)

Setcolor()

Drwline(G1,syasen[[1]])
Drwline(Gsall,Gsall2)
Drwline(Gyoko,Gtate)

#Drwline(syasen[[2]],0.5)

#九九の数値を表内に配置
Fontsize("ss")
for ( i in 1:9){
  for( j in 1:9){
     Putcell(hyou,i+1,j+1,"c",i*j)
  }
}

#表の1行と1列のセルに数字(1から9)を配置
Putrow(hyou,1,"d","","1","2","3","4","5","6","7","8","9")
PutcoL(hyou,1,"d","","1","2","3","4","5","6","7","8","9")

#左上隅セル内の文字配置
Fontsize("n")
Expr(c(Pt11[1],Pt11[2]),"n-2.5e","j")
Expr(c(Pt11[1],Pt11[2]),"s-1w0.5","i")

#表の上部中央にj列と表示
Letter(hyou_yoko_center[[1]]+c(0,2*hyou_yoko_center[[3]]+1),"n-2","$j$列")

#表の左中央にi行と表示
Letter(hyou_tate_center[[1]]-c(2*hyou_tate_center[[2]],0),"w-2","$i$行")

Endpicture(0)
Closefile()
####

■作表コマンド群は次を参照してください.(R版KETpicマニュアルでは未整備)
・KETpic v3.5.0 コマンド一覧for Scilab
「referencesci3_5_0.pdf」をダウンロード
・参考資料(使用例として)
「sakuhyou.pdf」をダウンロード



2011年7月 5日 (火)

直方体の有角透視図  二つの消点を用いてかく

■有角透視図の学習
直方体の有角透視図を二つの消点を用いてかく(二消点法)
教科書 土木製図(実教出版) 60ページを参考

■課題2(60ページ)  視点の位置をかえて、直方体の有角透視図を描く.
課題1  「toushizu110705_a_b.pdf」をダウンロード
課題2  「toushizu110705_c_d.pdf」をダウンロード

■課題の透視図例(R版KETpicにより作成)
各自で作成した透視図と比べてみよう.
「toushizu_rei_110705.zip」をダウンロード
(PDFファイルを圧縮)

Yukaku_toushizu_a1 Yukaku_toushizu_c1

 Yukaku_toushizu_b1  Yukaku_toushizu_d1


2011年6月 4日 (土)

50kgNレール断面図 の模写(土木製図)

【課題】 50kgNレール断面図 S=1:1の模写
工業026 土木製図(実教出版,平成21年2月)p71

・円と直線に接する半径rの円の描き方
・異なる半径の円の接続の方法
・勾配の表し方
など学ぶ.
半径500mm,300mmの円弧は円弧チェック用紙を活用して描いてみよう.
円弧定規に刻まれた線は円弧長10mm間隔で,円中心に向いている.
各円弧定規の幅は20mm.

●作図結果の確認用紙【実寸印刷して利用】(2011/06/25追記)
作図した図面の完成度を確認しましょう.実寸図を用意しました.
・大きくずれている場合には書き直しです.
・手で作図することが大切,確認用紙を写すことはNGです.
・この図は,作図習得度を確かめるために利用します.
.書き方が理解できなかった部分はメモしておく.
▼作成例(110628版) R版KETpicにより作成
「rail_JIS_50kgN_110628.pdf」をダウンロード
実寸印刷はインクジェットプリンターが良いようです.
レーザープリンターでは拡大縮小の調整が必要となるようです.たとえば,101%など.
この図面は,学習用として利用するものです.

【覚書メモ】
・センターライン(CL)記号はCとLを重ねて書く.
・Texファイルの書き出し部のR版KETpicのソースコード例
Openfile("zu_rail_50kgN.tex")
  Beginpicture("1mm")
    Drwline(RaildataR,RaildataL,1)
    Drwline(CL_hasen1_p,CL_hasen1_m,Hcyusin_hasen1,1)
    Letter(c(-30,76),"n","{\\it C}")
    Letter(c(-29.2,76),"n","{\\scalebox{1.2}{\\it L}}")
  Endpicture(0)
Closefile()

【確認してみよう】教科書の図面の中立軸の寸法に関する質問がありました.
JIS E1101に図が掲載されているようです.
JIS(日本工業規格)
▼普通レールの形状(JIS 50kg N)  
以下のパンフレット(JFE スチール株式会社)の7ページを参考にしてください.きれいな図面での確認できます.
http://www.jfe-steel.co.jp/products/katakou/catalog/d1j-003.pdf

■円弧チェック用紙 (原寸印刷して利用)
「enko_jyoug_110604.pdf」をダウンロード (R=80mm,300mm,500mm円弧定規)
「enko_jyougi300_470_500mm.pdf」をダウンロード  (R=300mm,470mm,500mm円弧定規)

■円弧チェック用紙(円弧定規R=500)作成のコード例(R版KETpic)
「ketpic110601.RData」をダウンロード

(各変数名の役割についいて,説明を受けてから利用してください)

#R版KETpic
#円弧定規の作成のソースコード例
#円弧の長さ10mmごとに円弧の中心方向が刻まれ,必要な縮尺と半径の大きさを記入
#R=500mm 縮尺1:1の例
rm(list=ls(all=TRUE))#ワークスペース上の全てのオブジェクトを消去する
setwd("C:/rwork")
load("ketpic110601.RData")
load("tgpack110328.RData")
R<- 500#半径mm
L<- 150#弦長mm
ang<- acos(L/2/R)
angL<- pi/2-ang
arcL<- R*angL
Memori<- 10#10mm
Num_memori<- floor(arcL/Memori)
Haba<- 20#定規の幅mm
theta<- asin(L/2/R)
begin_kaku<- pi/2-theta
end_kaku<- pi/2
C1<- Circledata(c(0,0),R,"N=400","R=c(begin_kaku,end_kaku)")#外側の円弧
pt_theta<- Memori/R
O<- c(0,0)
Pt_beg<- R*c(cos(begin_kaku),sin(begin_kaku))
Pt_Top<- c(0,R)
migi_line<- Listplot(Pt_beg,c(R*cos(begin_kaku),0))
C2_tmp<-  Circledata(c(0,0),R-Haba,"N=400")#内側の円
kouten<- Intersectcrvs(C2_tmp,migi_line)[[1]]
en2_beg_kaku<- atan(kouten[2]/kouten[1])
C2<-  Circledata(c(0,0),R-Haba,"N=400","R=c(en2_beg_kaku,pi/2)")#内側の円弧
pt_theta<- Memori/R
Pt_all<- list()
for (i in 1:Num_memori){
  ptn<- (pi/2-i*pt_theta)
  Pt_all<- c(Pt_all,list(R*c(cos(ptn),sin(ptn))))
}
Line_all<- list()
for (i in 1:Num_memori){
  Tmp<- Listplot(Pt_all[[i]],c(0,0))
  Line_all<- Flattenlist(list(Line_all,Tmp))
}
kouten_small_EN<- Intersectcrvs(C2,Line_all)
Line_all_kouten_s_EN<- list()
for (i in 1:Num_memori){
  tmp<- Listplot(kouten_small_EN[[i]],Pt_all[[i]])
  Line_all_kouten_s_EN <- Flattenlist(list(Line_all_kouten_s_EN,tmp))
}
Migi_side<- Listplot(Pt_beg,kouten)
Suiheisen<- Listplot(c(0,kouten[2]),kouten)
Right_G<- c(list(C1,C2,Line_all_kouten_s_EN ,Migi_side,Suiheisen))
Right_G<- Flattenlist(Right_G)
Left_G<- Reflectdata( Right_G,c(c(0,0),c(0,1)) )
kado<- c(L/2,Pt_Top[2])
H1<- Listplot(Pt_Top,kado)
V1<- Listplot(Pt_beg,kado)
Right_H1_V1<- c(list(H1,V1))
Left_H1_V1<- Reflectdata( Right_H1_V1,c(c(0,0),c(0,1)) )
T_R_L_H1_V1<- c(list(Right_H1_V1,Left_H1_V1))
Center_line<- Listplot(Pt_Top,c(0,R-Haba))
Total_G<- c(list(Right_G,Left_G,Center_line))
Pt_l3<- c(0,kouten[2])

Setwindow(c(-L/2,L/2),c(R-2*L/2,R))
WindispT()
WindispT(Total_G,T_R_L_H1_V1)

#以下,図ファイル(zu_enko_jyougi_test)の作成
#作成した図ファイルをTEXで処理 \input{zu_enko_jyougi_test.tex}
Setwindow(c(-L/2,L/2),c(kouten[2],Pt_Top[2]))
Openfile("zu_enko_jyougi_test.tex")
  Beginpicture("1mm")
    Drwline(Total_G,1)
    Drwline(T_R_L_H1_V1,0.5)
    Letter(Pt_l3,"n","R=500")
    Letter(Pt_l3,"s","1:1")
  Endpicture(0)
Closefile()

2011年5月 9日 (月)

学習用三角スケールの目盛(PDFファイル)

■学習用三角スケールの目盛(15cm)PDFファイル(学習用に利用して下さい)
KETpicで作成

図学,製図,測量実習などの学習用として役立てて下さい.
以下の4を推奨.

1) 三角スケール目盛05(8ミリインチ×0.5の太さ)
目盛(1/100,1/200,1/250,1/300,1/500,1/600 m)
「trianglescale_05.pdf」をダウンロード

2) 三角スケール目盛03(8ミリインチ×0.3の太さ)
目盛(1/100,1/200,1/250,1/300,1/500,1/600 m)
「trianglescale_03.pdf」をダウンロード

3) 三角スケール目盛05(8ミリインチ×0.5の太さ)
目盛(1/100,1/150,1/200,1/250,1/300,1/400,1/500,1/600 m)
「trianglescale100_150_600_all_05.pdf」をダウンロード

4) 三角スケール目盛り02(8ミリインチ×0.2の太さ)(推奨)
目盛(1/100,1/200,1/300,1/400,1/500,1/600 m 三角スケール風)
A4サイズに3セット含まれます.
「sansuke110415a.pdf」をダウンロード

※一度,三角スケール目盛り作成プログラムを作成すると,少し修正することで,任意の目盛り目盛りが作れます.
また,縮小目盛り(70.7%縮小など),インチ目盛,大工さんの使う寸目盛りなども作成することもできます.

印刷方法(1) 【インクジェットプリンタ】
インクジェットプリンタでは原寸の印刷ができるようです.

※OHPフイルムなどの透明のものに印刷して使うことも考えられます.
・インクジェットプリンタ用OHPフイルムに直接印刷
・A4用紙に印刷したものを,複写機(モノクロ)用のOHPフィルムに複写

▼印刷方法(2)
プリンターにより異なるかもしれませんが,情報センターのEPSON LP7700の場合,以下のように設定して印刷すると,ほぼ15cm幅の原寸目盛の印刷ができると思います.
他のプリンターでは適宜調整する必要があると思います.

Adobe Acrobat Reader9でファイルを開く
●ファイルメニューの印刷を選択
ページの拡大/縮小  (なし)
自動回転と中央配置 (チェック外す)
PDFのページサイズに合わせて用紙を選択 
(チェック外す)
EPSON LP-7700 のプロパティ
・基本設定
印刷品質 (きれい)
・レイアウト
出力用紙(A4)
拡大/縮小 (101%)

2011年4月26日 (火)

製図例の写図の学習(テンプレートの作成)

製図例の写図等の学習に役立てる目的で,いくつかのテンプレートをKETpicで作成します.

■街路標準構造図 一般横断面図 S=1:50
(土木製図(7 実教 工業026),製図例1,実教出版)の作図から学ぶこと
▼1:50 直線勾配とは
ここでは,、水平長さ(x)と高さ(y)との比(y/x)を勾配と呼び,歩道部の雨水は,1/50の傾きにより,雨水ます(道路に敷設されているので観察してみよう)に流入することになります.
▼1:50 放物線勾配とは
1:50は,車道の横断勾配は車道端部と車道中央部を直線で結んだ直線勾配を指している.
車道の横断勾配は車の流れ方向と直角の方向の勾配意味する.
道路幅員(道路の幅)方向の横断形状は、路面排水の目的もあり,中央部は端部より高くなっている.
放物線勾配は,車道端部と車道中央部を放物線で結ぶことを示している.つまり車道の横断形状の全体は,上に凸の一つの放物線形状をしている.

■放物線勾配に関する基本式とそのテンプレートを作る
ここでは次の放物線を考えることにする.

Ax2  (1式)

下図の横断勾配の形状に対応させて(1式)の係数a,bを算出する.
x=0のとき,y=c ゆえに,(1)式のb=c
x=ℓのとき,y=0 ゆえに,
Acl2_2

放物線勾配に関する基本式は次式となる.

Oudanhobutu
y=c(1-(1/ℓ^2)*x^2)

Houbutu

【デザインデータブック,8-16 曲線縦距および曲線長
(日本橋梁建設協会)を活用させていただきました】

■課題
  1:50放物線勾配として,車道の横断勾配線を描く.幅員ℓ=7835mmとする.
【解】 c/ℓ=1/50 ゆえに c=156.7mm
       放物線勾配に関する基本式は次式となる.
       y=156.7*[1-(1.63×10^-8)*x^2]

■応用課題
製図例の縮尺はS=1:50である.課題の1:50放物線勾配のテンプレートを作成する.
(CADで描くことができるが,工夫しながら手書きする作業は大切な学習となります)
▼作画例
「houbutu110426.pdf」をダウンロード

【注意】放物線勾配の場合,道路端部での接線勾配は直線勾配の2倍となる.

▼ソースコード例(R版KETpic)
実寸(縮尺図面にも対応)の図面を作成のために,TEXで読み込むための図ファイルの書き出し(出力)が必要になります.この図ファイルをTEXに読み込みPDFファイルを作成します.
以下のソースコードは図ファイルを書き出すまでのものです.

############
rm(list=ls(all=TRUE))#ワークスペース上の全てのオブジェクトを消去する
setwd("C:/rwork")
load("ketpic101207.RData")
load("tgpack110328.RData")

slope<- 1/50 #放物線勾配1:50
l<- 7835 #単位mm(幅員/2)
c<- l*slope #156.7mmになる
Setwindow(c(0,l),c(0,c))
WindispT()
#右車線の横断面
G<- Plotdata("c*(1-1/l^2*x^2)","x=c(0,l)")
G1<- Listplot(c(0,0),c(0,c))
G2<- Listplot(c(0,0),c(l,0))
WindispT(G,G1,G2,new=TRUE)
#左車線の横断面
Gref<- Reflectdata(G,c(c(0,0),c(0,1)))
Greft<- Translatedata(Gref,l,0)
Gl1<- Listplot(c(l,0),c(l,c))
Gl2<- Listplot(c(0,0),c(l,0))
WindispT(Greft,Gl1,Gl2,new=TRUE,color="blue")

#以下はTEXで読み込むための図ファイルの書き出し(出力)
#右車線の横断面の図ファイル出力
Openfile("zu_houbutusen1.tex") #図ファイル
  Beginpicture("0.02mm")#縮尺1:50に設定
    Drwline(G,G1,G2,0.2)#線の太さ0.2を設定
  Endpicture(0)
Closefile()

#左車線の横断面の図ファイル出力
Openfile("zu_houbutusen2.tex") #図ファイル
  Beginpicture("0.02mm")##縮尺1:50に設定
    Drwline(Greft,Gl1,Gl2,0.2)#線の太さ0.2を設定
  Endpicture(0)
Closefile()
############

■等角図だ円用テンプレートの作成(学習用として活用します)
R版KETpicで作成.
だ円の斜径の1.2247倍がだ円の長径になる   (sqrt(3)/sqrt(2)=1.2247)
「daen_template110425.pdf」をダウンロード

▼PDFファイルを実寸で印刷
・インクジェットプリンターでは実寸が得られました.
・現在使用中のモノクロレーザープリンターでは101%の拡大で印刷
・プリンターの機種により調整して印刷する.

2011年3月 8日 (火)

KETpicの本が出版されました(Scilab版)

CASTEX応用研究会編,KETpicで楽々TEXグラフ,イーテキスト研究所,2011.3.4

数学の教科書にでてくる図形が多数載っています.
工学系における作図にも利用できます.
プログラミングに慣れていない方でも,“描画してみよう”という気持ちがあれば,
少しずつ使えるようになると思います.
載っているサンプルで練習してみましょう.コマンドの機能は
使ってみると分かってくると思います.
図形に触れることの多い方々(学生さんの学習にも)に使って頂きたいソフトの一つです.

2011年3月 5日 (土)

KETpic 入門セミナーのご案内 平成 23 年 3 月 29 日(火)

国立木更津工業高等専門学校  KETpic 入門セミナー
R, Scilab による図入り教材作成ツール入門

  KETpic は,Scilab, R などの数式処理システムを用いて,
TEX 文書へ図を挿入することが簡単にできるパッケージです.
KETpic を用いれば,正確で美しい図の入った文書や教材を
思い通りに作成できます.本セミナーでは,Scilab 版または
R 版 KETpic を用いた教材や文書作りを体験できます.

1. 日時:平成 23 年 3 月 29 日(火)13:00~17:00
2. 会場:木更津工業高等専門学校  ネットワーク情報センター  演習室
3. 主催:木更津工業高等専門学校  基礎学系数学科
4. 共催:CASTeX 応用研究会(キャステック)
5. 日程:
  13:00 - 13:30    KETpicとは               高遠節夫(東邦大学)
  13:30 - 14:20    平面曲線の描画        金子真隆(木更津高専)
  14:30 - 15:20    空間曲線の描画        山下 哲(木更津高専)
  15:40 - 16:20    R 版 KETpicの紹介    高橋克夫(木更津高専)
  16:20 - 17:00    作図演習
6. 参加費:無料
7. 申し込み先:詳しくはここをご覧ください

深井戸(被圧地下水)の流線網を描く

■深井戸(被圧地下水)の流線網を描く
深井戸(被圧地下水)の流線網を描く例題をKETpicを用いて描く.
FAIR/GEYER/OKUN,Water and Wastewater Engineering Volume1. John Wiley &Sons 9-18ページの例題Example 9-3.

Flow_net_3

Flow_netd

■描画結果

・座標値の単位はメートル.
・横軸x,縦軸z
・井戸は座標(0,0)に置かれている.
・流線網から,赤色に内側の流線は井戸に流入する.
・P点はよどみ点(stagnation point)となり,その座標はP( -5.433168m, 0)と算出される.
・縦軸z>0 はベクトル場(矢印)を重ねて表示.
・縦軸z<0 は鉛筆で流線を作図する手法例.

■計算及び描画条件
●計算パラメータの値
・長さの単位ft(feet)
K<- 3.28*10^-4  fps(feet per second) 透水係数
L<- 140 ft 井戸から上流方向への距離
p0<- 30 ft 井戸での低下水頭
s<- p0/L  ピエゾ水頭の勾配 (井戸に現れる水面の高さをピエゾ水頭と呼ぶ)
m<- 40 ft 帯水層厚さ
w<- 56 ft  計算対象とする帯水層幅の1/2(幅は2wとする)
a<- 2*w*m ft^2  計算対象とする流下方向の「帯水層の断面積」
Q<- K*s*a ft^3/s 井戸の揚水量
q<- Q/m ft^3/ft/s 単位帯水層厚当たりの揚水量
v0<- K*s  ft/s 浸透層内の流速成分Darcyの法則v0=Ksにより算定.井戸のない条件でのDarcy速度(平行流)

・長さの単位m
unit<- 30.48/100 換算係数
p0<- 30*unit
L<- 140*unit
K<- 3.28*(10^-4)*unit
s<- p0/L
w<- 56*unit
m<- 40*unit
a<- 2*w*m
Q<- K*s*a
q<- Q/m
v0<- K*s

■上記パラメータを用いて算出したQ,q,v0の値
・井戸の揚水量
Q=0.31488f^3/s
・単位帯水層当たりの揚水量
q=0.007872f^3/f/s
・Darcy速度(平行流)
v0=7.028571e-05f/s

Q=0.008916409m^3/s
q=0.0007313327m^3/m/s
v0=2.142309e-05m/s

■流線描画の手順
●水理式
x方向への平行流の速度ポテンシャル
φ0(x,z)=v0x
井戸の吸い込みによる速度ポテンシャル
φr(x,z)=(q/2π)ln(R/r)=(q/4π)ln[(x2+z2)/r2]
二つを組み合わせた速度ポテンシャル
φ=φ0r =v0x+(q/4π)ln[(x2+z2)/r2]

ゆえに,以下の連立方程式より流線を求めることになる.
vx=∂φ/∂x=v0+(q/2π)[x/(x2+z2)]
vz=∂φ/∂z=qz/(2π)/(x2+z2)

●流線のプロットデータPDのDL
プロットデータをDLして描画してみよう.
ベクトル場などの矢印の配置方法を練習する.

■ソースコード例(R版KETpic)
1) tgpack110321.RDataの入手
「tgpack110321.RData」をダウンロード

2)流線データの入手
「idoryusen.zip」をダウンロード

####
rm(list=ls(all=TRUE))#ワークスペース上の全てのオブジェクトを消去する
setwd("C:/rwork")
load("ketpic101207.RData")
load("tgpack110321.RData")
#関数
Vx_Vz<- function(xz){#速度ベクトルの算出
  x<- xz[1]
  z<- xz[2]
  vx<- -v0-q*x/(2*pi)/(x^2+z^2)
  vz<- -q*z/(2*pi)/(x^2+z^2)
  return(list(vx,vz))
}

FlowlineArrow<- function(Str,size,color="black",width=1){
#流線上の矢印の表示
  for (i in 1:nrow(Str)){#i=1でNGになることの対応
    if (i==1){
      x<- Vx_Vz(Str[i,])[[1]]
      y<- Vx_Vz(Str[i,])[[2]]
      ArrowheadT(Str[i,],c(x,y),size,color=color,width=width)
    }else{
      ArrowheadT(Str[i,],Str,size,color=color,width=width)
    }
  }
}
#######
unit<- 30.48/100
p0<- 30*unit
L<- 140*unit
K<- 3.28*(10^-4)*unit
s<- p0/L
w<- 56*unit
m<- 40*unit
a<- 2*w*m
Q<- K*s*a
q<- Q/m
v0<- K*s

yodomiten<- -q/(2*pi*v0)
cat("P(x0)= ",yodomiten,"\n")

Setwindow(c(-30,55),c(-26,26))
WindispT()

#流線の描画
Ga1<- as.matrix(Readtextdata("flow_line_a1.txt","Cna=TRUE"))#流線データの読み込み
Ga2<- as.matrix(Readtextdata("flow_line_a2.txt","Cna=TRUE"))
Ga3<- as.matrix(Readtextdata("flow_line_a3.txt","Cna=TRUE"))
Ga4<- as.matrix(Readtextdata("flow_line_a4.txt","Cna=TRUE"))
Gb1<- as.matrix(Readtextdata("flow_line_b1.txt","Cna=TRUE"))
Gb2<- as.matrix(Readtextdata("flow_line_b2.txt","Cna=TRUE"))
Gb3<- as.matrix(Readtextdata("flow_line_b3.txt","Cna=TRUE"))
Gb4<- as.matrix(Readtextdata("flow_line_b4.txt","Cna=TRUE"))
Gb5<- as.matrix(Readtextdata("flow_line_b5.txt","Cna=TRUE"))
Gb6<- as.matrix(Readtextdata("flow_line_b6.txt","Cna=TRUE"))
Gb7<- as.matrix(Readtextdata("flow_line_b7.txt","Cna=TRUE"))
Gb8<- as.matrix(Readtextdata("flow_line_b8.txt","Cna=TRUE"))
Gb9<- as.matrix(Readtextdata("flow_line_b9.txt","Cna=TRUE"))#赤色流線

Ga<- Joingraphics(Ga1,Ga2,Ga3,Ga4)#流線
Gb<- Joingraphics(Gb1,Gb2,Gb3,Gb4,Gb5,Gb6,Gb7,Gb8)#赤色内の流線

Gb9_xminRow<- which(Gb9==min(Gb9[,1]))#赤色流線の前処理
GP9part<- Partcrv(Gb9[Gb9_xminRow,],Ptend(Gb9),Gb9)

WindispT(GP9part,new=TRUE,color="red",width=2)

#下部の放射線の描画
G1<- list()
for (i in 1:7){
  Tmp<- AngLine(c(0,0),-22.5*i,70)
  G1<- Flattenlist(list(G1,Tmp))
}
WindispT(Ga,Gb,G1,new=TRUE)

#下部の破線の描画
kankaku<- 2
k_line<- list()
ky<- -kankaku
step<- ky
while  (ky >= -26){
  k_line<- Flattenlist( list(k_line,Listplot(c(-50,ky),c(60,ky))) )
  ky<- ky+step
}
Tmp<- as.matrix(k_line)
for (i in 1 :(length(k_line))){
  K_line<- DashlineT(Tmp[i,],7)
  WindispT(K_line,new=TRUE)
}

rad_a<- Intersectcrvs(G1[[1]],k_line)#放射線と横平行線群との交点
rad_b<- Intersectcrvs(G1[[2]],k_line)
rad_c<- Intersectcrvs(G1[[3]],k_line)
rad_d<- Intersectcrvs(G1[[4]],k_line)
rad_e<- Intersectcrvs(G1[[5]],k_line)
rad_f<- Intersectcrvs(G1[[6]],k_line)
rad_g<- Intersectcrvs(G1[[7]],k_line)

Str_1<- Listplot(rad_a[[1]])#下部の流線
FlowlineArrow(Str_1,5,color="blue")#下部の流線上の矢印
Str_2<- Listplot(rad_a[[2]],rad_b[[1]])
FlowlineArrow(Str_2,5,color="blue")
Str_3<- Listplot(rad_a[[3]],rad_b[[2]],rad_c[[1]])
FlowlineArrow(Str_3,5,color="blue")
Str_4<- Listplot(rad_a[[4]],rad_b[[3]],rad_c[[2]],rad_d[[1]])
FlowlineArrow(Str_4,5,color="blue")
Str_5<- Listplot(rad_a[[5]],rad_b[[4]],rad_c[[3]],rad_d[[2]],rad_e[[1]])
FlowlineArrow(Str_5,5,color="blue")
Str_6<- Listplot(rad_a[[6]],rad_b[[5]],rad_c[[4]],rad_d[[3]],rad_e[[2]],rad_f[[1]])
FlowlineArrow(Str_6,5,color="blue")
Str_7<- Listplot(rad_a[[7]],rad_b[[6]],rad_c[[5]],rad_d[[4]],rad_e[[3]],rad_f[[2]],rad_g[[1]])
FlowlineArrow(Str_7,7,color="red",width=1)
Str_8<- Listplot(rad_a[[8]],rad_b[[7]],rad_c[[6]],rad_d[[5]],rad_e[[4]],rad_f[[3]],rad_g[[2]])
FlowlineArrow(Str_8,6,color="blue")
Str_9<- Listplot(rad_a[[9]],rad_b[[8]],rad_c[[7]],rad_d[[6]],rad_e[[5]],rad_f[[4]],rad_g[[3]])
FlowlineArrow(Str_9,6,color="blue")
Str_10<- Listplot(rad_a[[10]],rad_b[[9]],rad_c[[8]],rad_d[[7]],rad_e[[6]],rad_f[[5]],rad_g[[4]])
FlowlineArrow(Str_10,6,color="blue")
Str_11<- Listplot(rad_a[[11]],rad_b[[10]],rad_c[[9]],rad_d[[8]],rad_e[[7]],rad_f[[6]],rad_g[[5]])
FlowlineArrow(Str_11,6,color="blue")
Str_12<- Listplot(rad_b[[11]],rad_c[[10]],rad_d[[9]],rad_e[[8]],rad_f[[7]],rad_g[[6]])
FlowlineArrow(Str_12,6,color="blue")
Str_13<- Listplot(rad_c[[11]],rad_d[[10]],rad_e[[9]],rad_f[[8]],rad_g[[7]])
FlowlineArrow(Str_13,6,color="blue")

Ppoint<- AngLine(Str_7[7,],90,4)#90度と4mを指定
Pten<- Ppoint[2,]-c(0,2)
WindispT(Ppoint,Pten,pch=19,new=TRUE,color="red",width=2)#よどみ点(P)pch=19黒丸.pch=21白丸

GStr<- Joingraphics(Str_2,Str_3,Str_4,Str_5,Str_6,Str_7,Str_8,Str_9,Str_10,Str_11,Str_12,Str_13)

WindispT(GStr,new=TRUE)
WindispT(Str_7,new=TRUE,width=2)
WindispT(Str_7,color="red",new=TRUE)#下部の赤線

#速度ベクトル場(矢印)の表示
xmesh<- c(seq(-32,52,by=4))
zmesh<- c(seq(4,26,by=4))

for (j in 1:length(zmesh)){#z axis
  for (i in 1:length(xmesh)){#x axis
   x<- xmesh[i]
   z<- zmesh[j]
   Tmp<- Vx_Vz(c(x,z))
   vx<- Tmp[[1]]
   vz<- Tmp[[2]]
   ArrowlineT(c(x,z),c(x,z)+100000*c(vx,vz),4,color="blue","tl")
  }
}
Puttext(c(53,-28),"11.03.03 K.takahashi",pos=2,cex=0.8)
Puttext(Pten,"P",pos=2,cex=1)

2011年3月 4日 (金)

一様水路の不等流(開水路)

■課題 一様水路の不等流(開水路)に関する描画
一様水路(prismatic channel):水路断面が不変で一定の底こう配をもった水路.

■用語の概要
参考図書(主なものを記載.古い図書が多くなりましが,他に新しい参考図書は多くあると思います)
1) 日野幹雄著 明解 水理学 ,平成6年10月,丸善(多くを参照)
2) 大西外明著 最新 水理学Ⅰ,Ⅱ 1982.10 ,森北出版
3) 椿東一郎著 水理学Ⅰ,Ⅱ 1973.2 ,森北出版
4) 椿東一郎・荒木正夫共著 水理学演習上,下,1984.2 ,森北出版
5) 岩佐義朗著  水理学,昭和46 年2月,朝倉書店
6) 岩崎敏夫著 水理学例題演習 ,昭和43 年4月,コロナ社
7) 岩崎敏夫著 応用水理学例題演習,昭和56 年6月,コロナ社
8) 本間仁著     標準水理学,昭和44 年3月,丸善
9) 本間仁著   水理学,昭和42 年5月,丸善
10) 鈴木幸一著 水理学演習2005.2,森北出版
11) VEN TE CHOW,OPEN-CHANNEL HYDRAULICS,McCRAW-HILL KOGAKUSHA LTD

1 定常流(steady flow):流速,水深,流水断面積等が時間的に変化しない流れ.
・等流(uniform flow):流れ方向に流速,水深等が変化せず,水面勾配と水路勾配が同じ流れ
・不等流(non-uniform flow,varied flow):流速,水深等が流れ方向に変化する流れ.

2 非定常流(unsteady flow):流速,水深,流水断面積等が時間的に変化する流れ.
・等流(unsteady uniform flow):まれ.
・不等流(unsteady varied flow):水面上を伝播する波がある流れ.

▼開水路における流れの分類(Ven Te Chow,Open-Chanel Hydraulics, p.7 より)
A. Steady flow
.  1. Uniform flow
  2. Varied Flow
      a. Gradually varied flow
      b. Rapidly varied flow
B. Unsteady flow
  1. Unsteady uniform flow(rare)
  2. Unsteady flow(i.e., unsteady varied flow)
      a. Gradually varied unsteady flow
      b. Rapidly varied unsteady flow

3 比エネルギー
 水路床から測ったその断面での単位体積重量あたりのエネルギーを水頭(水柱高さ換算)で表したものを比エネルギー(specific energy)という.水路底が水平及び水路底の勾配が小さい緩こう配水路の場合の比エネルギー(E)は,
   E=v2/2g + h
となる.h:水深,v:断面平均流速,g:重力加速度
幅の広い長方形断面水路の単位幅あたりの流量をqとすると,流速vは
     v=q/h
したがって,単位幅流量qが与えられた場合の比エネルギーは

  E(h,q)=q^2/(2g)・1/(h^2) + h

この関係を示す図を比エネルギー図という.

4 限界水深(critical depth),限界流速(critical velocity)
  一定の単位幅流量qに対し,比エネルギーが最小となる(比エネルギーがこれ以下にはなり得ない限界状態の)水深hcを限界水深とよぶ.

  Hc

この水深で比エネルギーが極小値になる.この場合は最小値でもある.

  限界状態での速度を限界速度vcとよぶ.
vc=√ghc    つまり   vc=sqrt(ghc)

・限界水深ではFroude数(フルード数)(Fr)は1である.
・限界水深は比エネルギーの2/3である.
・一定単位幅流量qに対する最小比エネルギーEc(q)は,以下の式で表される.

  Ecq

●流量をある一定値に保った場合,一つの比エネルギーEに対し,流れは限界水深の状態を除いて,常に二つの水深h1,h2をとることが可能である.一つは限界水深より大きく,他は限界水深より小さい水深.
   h1>hc>h2 さい.
  水深h1:常流→深く遅い流れ(Fr<1)
 水深h2:射流→浅く速い流れ(Fr>1)

5 流量図(discharge diagram)
  E(h,q)=q^2/(2g)・1/(h^2) + h より
  q^2=2g(E-h)h^2
これは,比エネルギーEをある一定値に保ったときの単位幅流量qと水深hの関係(流量図) 
Eが一定のもとで単位幅流量を最大qcとする水深は,
   Hce
となる.

6 等流水深(normal depth)と限界水深(critical depth) ,限界勾配(critical slope)
 等流水深:一様断面の水路が十分長く続いていると,水深も一様な状態に達する.この状態を等流状態とよび,このような流れを等流(uniform flow)という.
このような平衡状態に達した水深を等流水深(normal depth)h0という.Manningの平均流速公式から,

  H0

  i0:水路勾配
したがって,ある単位幅流量qに対する等流水深h0は水路勾配i0が急になるにつれて減少する.この関係を次の記号で示す.
    i0  ↑ :  h0  ↓ 

一方,水路勾配i0の増加とともにFr は増大し,ある勾配で Fr は1となる.
    i0  ↑ :  Fr  ↑ 

 ある単位幅流量qに対する等流がちょうどFroude数1となる水路勾配,あるいは等流水深が限界水深となる(h0=h)水路勾配を限界勾配(critical slope)ic という.
  Ic

  Fr =1となる限界水深hは,次式から,水路勾配  i0  に無関係に,単位幅流量qのみにより一義的に定まることが示される.

  Hcq2

▼課題  比エネルギー図(specific energy diagram)<h~E図>の概略図を描く.
参考図書1)  p.109

Sediagram_2

■ソースコード例

rm(list=ls(all=TRUE))#ワークスペース上の全てのオブジェクトを消去する
setwd("C:/rwork")
load("ketpic101207.RData")
load("tgpack110321.RData")

Setwindow(c(0,16),c(0,16))
qval<- c(50,100,150,200,250)#qの値の設定
hs<- c(0.30,0.6,0.9,1.2,1.5)#表示開始のhの設定
he<- c(8.8,8.7,8.6,8.5,8.4)#表示終了のhの設定

Zu_E_h<- list()# E縦軸_h横軸 基本式用
Zu_h_E<- list()# h縦軸_E横軸 描画用

#h-E図のPD
for (i in 1:5){
  Assignset("q",qval[i],"g",980,"hs",hs[i],"he",he[i])
  G<- Plotdata(Assign("q^2/(2*g)*(1/h^2)+h"),Assign("h=c(hs,he)"))#E-h図基本式用
  G1<- Reflectdata(G,c(c(0,0),c(1,1)))#h-E図 線対称移動関数の利用
  Zu_h_E<- Flattenlist(Zu_h_E,G1)#list構造の平準化
  Zu_E_h<- Flattenlist(Zu_E_h,G)
}

G45<- Listplot(c(0,0),c(10,10))#45度の線
WindispT(Zu_h_E,G45)#h-E図の描画.45度の線の描画

Pt1<- Ptcrv(2,Zu_h_E[[1]])#曲線PDのn番目の節点を取得
Pt2<- Ptcrv(4,Zu_h_E[[2]])
Pt3<- Ptcrv(7,Zu_h_E[[3]])
Pt4<- Ptcrv(12,Zu_h_E[[4]])
Pt5<- Ptcrv(20,Zu_h_E[[5]])
#矢印付き引出し線の描画(関数LeaderLineArの利用)
LL1<- LeaderLineAr(Pt1,c(-1,-1),2.5,arsize=1.5,color="black",width=1,"tf")
LL2<- LeaderLineAr(Pt2,c(-1,-1),2.5,arsize=1.5,color="black",width=1,"tf")
LL3<- LeaderLineAr(Pt3,c(-1,-1),2.5,arsize=1.5,color="black",width=1,"tf")
LL4<- LeaderLineAr(Pt4,c(-1,-1),2.5,arsize=1.5,color="black",width=1,"tf")
LL5<- LeaderLineAr(Pt5,c(-1,-1),2.5,arsize=1.5,color="black",width=1,"tf")

Puttext(LL1[2,],"50",pos=4)
Puttext(LL2[2,],"100",pos=4)
Puttext(LL3[2,],"150",pos=4)
Puttext(LL4[2,],"200",pos=4)
Puttext(LL5[2,],"q=250",pos=4)
Puttext(c(3,10),expression(paste("E =", frac(q^2,2*g),"・",frac(1,h^2),"+h",sep="")))#数式表示

Puttext(c(15.5,-0.2),"E",pos=3)#横軸ラベル
Puttext(c(0,11),"h",pos=2)#縦軸ラベル

#角マークの描画
kaku<- Anglemark(c(1,0), c(0,0),c(1,1),2)#角マークのPD
WindispT(kaku,new=TRUE)#角マークの表示
ArrowheadT(Ptend(kaku),kaku,size=1.5,color="black",width=1)

cyuten<- Bowmiddle(kaku)#円弧の中心座標を取得
Puttext(cyuten,"45°",pos=4)#円弧の中心に45°を表示


■課題 一様水路の不等流(開水路)に関する描画 【参考図書6) p.126 [問題1]】
  広矩形断面で一様な水路勾配i0=1/1320 ,B=200m,Q=1000m^3/sec,
h=6.5m,n=0.023として,背水曲線を描く.KETpicの微分方程式の解を求める関数(Deqplot)を用いる.(Bresseの背水公式によらない)
手順
1)  等流水深h0 限界水深hcを算出する.

  H0hc

ゆえに, h0 =2.358m  hc=1.410mと算出される.

2)  下記の不等流に関する基本方程式を解く.

    Dhdx

■水面形の描画(R 版KETpicによる)

Suimenkeim1

■ソースコード例
//Scilab版
clear;
Ketlib=lib('C:/sciwork/ketpicsciL5');
Ketinit();
Setwindow([-10000,0],[0,15]);
Setscaling(500);
cd('C:/sciwork');
g=9.8;
i0=1/1320;
h1=6.50;
B=200;
Q=1000;
n=0.023;
h0=((Q^2)*(n^2)/((B^2)*i0))^(3/10);
hc=((1.1*Q^2)/(g*(B^2)))^(1/3);
Assignset('h0',h0,'hc',hc);
Eq='h`=i0*(1-(h0/h)^3)/(1-(hc/h)^3)';
G=Deqplot(Assign(Eq),'x=[0,-10000]',0,h1);//水深のPD
Gi=-i0*G(:,1);//基準面から水路床までの高さ(z).水路勾配*L
Gislope=[G(:,1) Gi];//列ベクトルの結合.zのPD
G2=G(:,2)+Gi;//列ベクトルの和  基準面から水面までの高さ(z+h)
GT= [G(:,1) G2(:,1)];//基準面から水面までの高さ(z+h)のPD
Windisp(G,GT,Gislope);//水深(h), z+h  ,  zの描画

//3種のPDデータのファイル出力. R版KETpicで読み込み,描画利用する.
Writetextdata('futouG.txt',G);//カレントディレクトリの出力
Writetextdata('futouGT.txt',GT);
Writetextdata('futouGislope.txt',Gislope);

■ソースコード例(R版)
1)PDデータの入手
「suimendata.zip」をダウンロード

rm(list=ls(all=TRUE))#ワークスペース上の全てのオブジェクトを消去する
setwd("C:/rwork")
load("ketpic101207.RData")
load("tgpack110321.RData")

Setwindow(c(-10000,0),c(0,12))
sc<- 600
Setscaling(sc)
WindispT()
G1<- as.matrix(Readtextdata("futouG.txt","Cna=TRUE"))
G2<- as.matrix(Readtextdata("futouGislope.txt","Cna=TRUE"))
G3<- as.matrix(Readtextdata("futouGT.txt","Cna=TRUE"))
Tmp<- G2[,2]-matrix(0.5,nrow=nrow(G2))
Bs<- cbind(G2[,1],Tmp)
G2start<- Ptstart(G2)
G2end<- Ptend(G2)
Bsstart<- Ptstart(Bs)
Bsend<- Ptend(Bs)
L1<- Listplot(G2start,Bsstart)
L2<- Listplot(G2end,Bsend)
Soko<- Enclosing(list(G2,L2,Invert(Bs),Invert(L1)))#水路底のPD

WindispT(G3,color="blue",width=4,new=TRUE)#水面形
WindispT(list(col="chocolate1",border="white",density=20,Soko),new=TRUE)
WindispT(G2,color="chocolate1",width=3,new=TRUE)
Puttext(c(-1500,11.5*sc),"広矩形断面での一様な水路勾配における背水曲線",pos=2)
Puttext(c(-10,7*sc),"h=6.5m",pos=2)
Tmp<- Ptend(G1)
Puttext(c(-10000,10.5*sc),paste("h=",round(Tmp[2],digits = 3),"m"),pos=4)
Puttext(c(-10,0*sc),"単位m",pos=2)
Puttext(c(-3000,2*sc),"水路勾配(i0)=1/1320, B=200m, Q=1000m^3/sec",pos=2)
Puttext(c(-7000,1*sc),"h1=6.5m, n=0.023",pos=2)
Puttext(c(-10,-0.5*sc),"11.03.21 K.takahashi",pos=2,cex=0.8)

より以前の記事一覧

無料ブログはココログ
2017年8月
    1 2 3 4 5
6 7 8 9 10 11 12
13 14 15 16 17 18 19
20 21 22 23 24 25 26
27 28 29 30 31