Created
October 5, 2016 21:28
-
-
Save dwoll/a3a0e70edaa4d5b7abce7ad55342ddb9 to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
## DJV 1 target, see http://imgur.com/LHiZFeG | |
## unit x y -> cm | |
## points on outer polygon (ring 1) | |
djv1 <- read.table(file=textConnection("x y | |
25.4 67.9 | |
21.2 59.3 | |
27.5 49.4 | |
36.7 46.2 | |
50.4 46.2 | |
57.2 47.3 | |
58.9 55.7 | |
56.6 65.1 | |
40.2 65.7 | |
30.9 66.3"), header=TRUE) | |
## assumed shot distribution | |
## bivariate normal centered at ring 10 + bias c(2, 3) cm | |
## about sdX=1.5 MOA sdY=2 MOA -> sdX=9 and sdY=12 cm | |
## correlation = 0.2 | |
ctr10 <- c(40.2, 55.7) # target center coords | |
mu <- ctr10 + c(2, 3) # systematic accuracy bias | |
sdX <- 9 | |
sdY <- 12 | |
corXY <- 0.2 | |
covXY <- sdX*sdY*corXY | |
sigma <- cbind(c(sdX^2, covXY), c(covXY, sdY^2)) # covariance matrix | |
## define polygonal integration region covering DJV 1 ring 1 | |
## install.packages(c("shotGroups", "polyCub")) | |
## install.packages("gpclib", type="source") | |
library(polyCub) # for polyCub.exact.Gauss() | |
library(gpclib) # for as.gpc.poly() | |
djv1Pts <- data.matrix(subset(djv1, select=c(x, y))) # just x, y coords | |
djv1Poly <- as(djv1Pts[chull(djv1Pts), ], "gpc.poly") # convert to gpc polygon | |
gpclibPermit() # necessary for running polyCub.exact.Gauss() | |
## integrate shot distribution over polygon = DVJ 1 ring 1 | |
polyCub.exact.Gauss(djv1Poly, mean=mu, Sigma=sigma, plot=FALSE) | |
# 0.5431248 | |
## approximate DJV 1 ring 1 polygon with elliptical integration region | |
S <- cbind(c(18^2, 0), c(0, 10^2)) | |
e <- solve(S) | |
## integrate shot distribution over ellipse | |
library(shotGroups) # for pmvnEll(), drawEllipse() | |
pmvnEll(r=1, sigma=sigma, mu=mu, e=e, x0=ctr10) | |
# 0.4977843 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment