#############################
#
#  レイヤーモデル
#
#############################
#
#  アルゴリズム
#
#  ray の各ステップごとの変数
#          仰角(in radian)
#          レイヤ( 0.1 m = 10cm = 1dm ごとを想定。)
#          位置(x)
#      これらをリストにしておく
#      object <- list( angle=0.1, layer=10, pos=0 )
#
# 屈折を光線が境界面に達するごとに計算する
#     ただし、全反射になったら上下の向きを変える。
#
#############################
nextstep <- function( obj, depth=0.1 ){
#
  if( obj$angle > 0 ){
    upward    <- T
    nextlayer <- obj$layer + 1
  }else{
    upward    <- F
    nextlayer <- obj$layer - 1
  }
#
  fracin    <- frac(obj$layer)
  fracout   <- frac(nextlayer)
  outsin    <- cos( obj$angle ) * fracin / fracout 

  if( outsin > 1 ){
    return( list( angle=     -obj$angle, layer=obj$layer, pos=obj$pos + depth * tan( pi/2 - abs(obj$angle) )))
  }else{
    if( upward ){
      return( list( angle= acos(outsin), layer=nextlayer, pos=obj$pos + depth * tan( pi/2 - abs(obj$angle) )))
    }else{
      return( list( angle=-acos(outsin), layer=nextlayer, pos=obj$pos + depth * tan( pi/2 - abs(obj$angle) )))
    }
  }
}

###################################
#
# 上方蜃気楼
#

frac <- function(n){
  return( 1 + ( 0.00027 + 0.00006*cos( 7*pi/20 + n / 1000 * 2 * pi ) * exp( - n/1000) ) * exp( - n / 20000 ))
}

# png("Upper_Mirage.png",width=800,height=600)
par(fig=c(5,40,0,100)/100)
plot(frac(1:2000),1:2000, ylim=c(0,2000), type='l',
     xlab="Refraction Index", ylab="Height[dm]")

par(fig=c(30,100,0,100)/100, new=T)
plot(0,0,type='n', xlim=c(0,30000), ylim=c(0,2000),
     xlab="Distance[dm]", ylab="")

for( angle in seq(0.02,-0.02,-0.0005)){
  object <- list( angle=angle, layer=10, pos=0 )
  for( i in 1:2000 ){
    object <- nextstep(object)
    par(mfg=c(1,1))
    points(object$pos, object$layer, pch='.')
    if( object$layer < 0 ) break
  }
}
# dev.off()

###################################
#
# 下方蜃気楼
#

frac <- function(n){
  return( 1 + ( 0.00027 - 0.00010*cos( -2*pi/20 + n / 400 * 2 * pi ) * exp( - n/100) ) * exp( - n / 20000 ))
}

# png("Lower_Mirage.png",width=800,height=600)
par(fig=c(5,40,0,100)/100)
plot(frac(1:600),1:600, ylim=c(0,600), type='l',
     xlab="Refraction Index", ylab="Height[dm]")

par(fig=c(30,100,0,100)/100, new=T)
plot(0,0,type='n', xlim=c(0,5000), ylim=c(0,500),
     xlab="Distance[dm]", ylab="")

for( angle in seq(0.015,-0.015,-0.0005)){
  object <- list( angle=angle, layer=40, pos=0 )
  for( i in 1:600 ){
    object <- nextstep(object)
    points(object$pos, object$layer, pch='.')
    if( object$layer < 0 ) break
  }
}
# dev.off()