不管有多烂，总是需要用的，备份一下 wellbore trajectory computation

我确定用不了三天，我就会忘记了这些东西是怎么捣鼓出来的，特别是乱七八糟的赋值与命名，所以，一定需要备份一下这些个代码，以后需要用的时候直接来扒拉。。。

1.采用校正角平均方法计算井斜

书名=定向钻井设计与计算 （第二版）

作者=韩志勇编

出版日期=2007.11

出版社=中国石油大学出版社

具体算法：校正角平均法见P75页

要特别注意角度与弧度之间的换算

2.内插值（斜深）计算

采用内插一个斜深，得到井斜与方位，然后用上面的校正角平均法重新计算一次，为什么要用这么没道理的方法呢，因为书上讲的那些个方法我算出来的结果差别很大，又找不到原因。。。

3.最简单的输入垂深返回斜深

给一个垂深，按最简单的直线模型算一个斜深、井斜与方位。

<-----------------

以下是这堆烂代码

初始化数据，主要是要避免参与运算时，第一行数据中方位角不要为0

init_data <- function(data){
  if (data[1,1]==0){
    data <- data[-1,]
  }
  data <- rbind(c(0,0,data[1,3]),data)
  return(data)
}

井斜计算-校正角平均法

trajcal <- function(data,H=0,S=0,N=0,E=0){
  depth=data[,1]
  dev=data[,2]
  az=data[,3]
  tmp <- c(depth=0,dev=0,az=az[1],H=H,N=N,E=E,V=0,allaz=0)
  tmp1<- NULL
  for(i in 2:length(depth)){
  delta_len <- depth[i]-depth[i-1]
  delta_dev <- (dev[i]-dev[i-1])*pi/180
  mean_dev <- 0.5*(dev[i]+dev[i-1])*pi/180
  delta_az <- az[i]-az[i-1]
  if (delta_az>180){
      deltaaz <- (deltaaz-360)pi/180
      mean_az <- (az[i]+az[i-1]-360)pi/360
  }else if(deltaaz<(-180)){
      deltaaz <- (deltaaz+360)*pi/180
      meanaz <- (az[i]+az[i-1]+360)pi/360
  }else{
      deltaaz <- deltaazpi/180
      mean_az <- (az[i]+az[i-1])*pi/360
  }


delta_H <- (1-delta_dev^2/24)*delta_len*cos(mean_dev)


deltaS <- (1-deltadev^2/24)delta_lensin(mean_dev)

delta_N <- (1-(delta_dev^2+delta_az^2)/24)*delta_len*sin(mean_dev)*cos(mean_az)
delta_E <- (1-(delta_dev^2+delta_az^2)/24)*delta_len*sin(mean_dev)*sin(mean_az)

H <- H+delta_H
N <- N+delta_N
E <- E+delta_E
V <- sqrt(N^2+E^2)

if (N>0){
 allaz <- atan(E/N)*180/pi
  }else if(N<0){
    allaz <- atan(E/N)*180/pi+180
}
if (allaz<0){
  allaz <- allaz+360
}

tmp1 <- cbind(depth=depth[i],dev=dev[i],az=az[i],H,N,E,V,allaz)
tmp <- rbind(tmp,tmp1)  


}
  return(tmp)
}

内插,给一个或一堆井深，算出对应的垂深与坐标，这样就可以把一堆傻傻的数据换算到垂深尺度上了。

trajint <- function(data,indepth){
  #调用trajcal函数进行初始的井斜计算
  result_data <- trajcal(data)
  tmp <- NULL
  for (i in 1:length(indepth[, 1])) {
#找到插值运算的第一个井深     
      x <- which(data[,1]>indepth[i,1])[1]
      if (is.na(x)){
        break
      }else{           
      deltaL <- indepth[i,1]-data[x-1,1]
      alpha <- data[x-1,2]+deltaL*(data[x,2]-data[x-1,2])/(data[x,1]-data[x-1,1])
  deltaaz <- data[x,3]-data[x-1,3]
az <- data[x-1,3]+deltaL*deltaaz/(data[x,1]-data[x-1,1])

#按这个井斜方位计算一次   
  delta_alpha <- (alpha-data[x-1,2])*pi/180
  mean_alpha <- 0.5*(alpha+data[x-1,2])*pi/180
  delta_az <- az-data[x-1,3]
  if (delta_az > 180){
      delta_az <- (delta_az-360)*pi/180
      mean_az <- (az+data[x-1,3]-360)*pi/360
  }else if(delta_az < (-180)){
      delta_az <- (delta_az+360)*pi/180
      mean_az <- (az+data[x-1,3]+360)*pi/360
  }else{
      delta_az <- delta_az*pi/180
      mean_az <- (az+data[x-1,3])*pi/360
  }

    delta_D <- (1-delta_alpha^2/24)*deltaL*cos(mean_alpha)
#     delta_S <- (1-delta_dev^2/24)*delta_len*sin(mean_dev)
    delta_N <- (1-(delta_alpha^2+delta_az^2)/24)*deltaL*sin(mean_alpha)*cos(mean_az)
    delta_E <- (1-(delta_alpha^2+delta_az^2)/24)*deltaL*sin(mean_alpha)*sin(mean_az)
  
    allH <- result_data[x-1,4]+delta_D
    allN <- result_data[x-1,5]+delta_N
    allE <- result_data[x-1,6]+delta_E
    allV <- sqrt(allN^2+allE^2)

   if (allN>0){
     allAZ <- atan(allE/allN)*180/pi
      }else if(allN<0){
        allAZ <- atan(allE/allN)*180/pi+180
    }
    if (allAZ<0){
      allAZ <- allAZ+360
    }
    tmp <- rbind(tmp, cbind(depth=indepth[i,1],dev=alpha,az=az,TVD=allH,N=allN,E=allE,V=allV,allAZ=allAZ,indepth[i,]))
      }
 }
    return(tmp)
}

插入垂深，反算斜深与坐标

trajinttvd <- function(data,tvd){
  #调用trajcal函数进行初始的井斜计算
  result_data <- trajcal(data)
  x <- which(result_data[,4]>tvd)[1]
  deltaD <- tvd-result_data[x-1,4]
  D <- result_data[x,4]-result_data[x-1,4]
  depth <- result_data[x-1,1]+(result_data[x,1]-result_data[x-1,1])*delta_D/D
  alpha <- result_data[x-1,2]+(result_data[x,2]-result_data[x-1,2])*delta_D/D
  az <- result_data[x-1,3]+(result_data[x,3]-result_data[x-1,3])*delta_D/D
  result <- cbind(depth,alpha,az)
  return(result)
}

范例

setwd("~/code/R_code/hf")

setwd("E:/项目/")
awb <- read.table("23.txt")
awb[1,3] <- 0

data <- awb

colnames(awb)<-c("Depth","Dev","Az")

初始化数据

awb1 <- init_data(awb)

井斜计算

abcd <- trajcal(awb1)

利用已有的井斜数据，内插值（斜深）计算

xrf <- read.csv("日报表.csv")
xrf_data <- trajint(awb1,xrf)

输入垂深，返回该点的斜深

tvd2depth <- trajint_tvd(awb1,3500)