Print

Please login or register.
1 Hour 1 Day 1 Week 1 Month Forever
Login with username, password and session length

[Liral]

Update: I have nearly rewritten one R script file to satisfy the conditions I had stipulated but need a couple files I can't find on the thread to run it.  Once I can run the script, I can fix the errors I have surely introduced because I have learned R just to make these changes.  I have also rewritten the code to load almost all data from .csv files.

Rewritten

Code: Rewritten

#use computationally inexpensive functions to manage

shipnames <- c("glimmer","brawlerlp","vanguard","tempest","medusa",
               "hammerhead","enforcer","dominator","fulgent","brilliant",
               "radiant","onslaught","aurora","paragon","conquest",
               "champion")

#get a filename
filename <- function(x) { 
    paste("optimizeweaponsbytime", 
          shipnames[x], 
          "allweaponswithacc.txt", 
          sep = "") 
}


#read a file
readfile <- function(filename) { 
    read.csv(filename, sep = "\t") 
}


#methods to get data tables

getFormattedTTKTable <- function(table) {
    table <- cbind(table, (1 / (table$Time / mean(table$Time)) - 1))
    table[,3] <- sprintf("%0.1f%%", table[,3] * 100)
    table <- table[order(table$Time),]
    table[,2] <- round(table[,2], digits=1)
    return(table)
}


getTTKTable <- function(dataFrame, time, label) {
    table <- aggregate(time, data = dataFrame, FUN = mean, na.rm = TRUE)
    table <- getFormattedTTKTable(table)
    colnames(table) <- c(label, "Avg. time to kill", "TTK speed score")
    return(table)
}


getGunAndLargeMissileTable <- function(dataFrameSorted) {
    table <- aggregate(Time ~ Largemissiles + Largeguns + Mediumguns, 
                       data = dataFrameSorted, 
                       FUN = mean, 
                       na.rm = TRUE)
    table <- getFormattedTTKTable(table)
    colnames(table) <- c("Large missiles", 
                         "Large guns",
                         "Medium guns", 
                         "Avg. time to kill", 
                         "TTK speed score")
    return(table)
}


#method to write data tables to the screen
writeTable <- function(table, filename) {
    write.table(table, file = filename, row.names = FALSE, sep = "/t")
}


getDataFrame <- function(data) {
    dataFrame <- data.frame()
    
    for (i in 1:length(shipnames)) {
        shipnamecolumn <- matrix(data=shipnames[i], ncol=1, nrow=7938)
        mainDataFrame <- rbind(dataFrame, 
                               cbind(shipnamecolumn, readfile(filename(i))))
    }

    colnames(dataFrame) <- c("Ship", 
                             "Time", 
                             "Largemissile1",
                             "Largemissile2", 
                             "Mediummissile1",
                             "Mediummissile2",
                             "Largegun1", 
                             "Largegun2", 
                             "Mediumgun1",
                             "Mediumgun2")
    
    return(dataFrame)
}


getDataFrameSorted <- function(data, dataFrame) {
    dataFrameSorted <- data.frame()
    
    #order large missiles, medium missiles, and large guns alphabetically by row 
    #and merge cells
    for (i in 1:length(dataFrame[,1])){
        largeMissiles <- paste(t(apply(dataFrame[i,3:4], 1, sort))[2],
                               t(apply(dataFrame[i,3:4], 1, sort))[1], 
                               sep = " ")
        largeGuns <- paste(t(apply(dataFrame[i,7:8], 1, sort))[1],
                           t(apply(dataFrame[i,7:8], 1, sort))[2], sep = " ")
        mediumGuns <- paste(t(apply(dataFrame[i,9:10], 1, sort))[1],
                            t(apply(dataFrame[i,9:10], 1, sort))[2], sep = " ")
        toBind <- c(dataFrame[i,1], dataFrame[i,2], largeMissiles, largeGuns, 
                    mediumGuns)
        dataFrameSorted <- rbind(dataFrameSorted, toBind)
    }
    
    colnames(dataFrameSorted) <- c("Ship", 
                                   "Time", 
                                   "Largemissiles", 
                                   "Largeguns", 
                                   "Mediumguns")
                                       
    dataFrameSorted$Time <- as.double(dataFrameSorted$Time)
    
    return(dataFrameSorted)
}


#analysis with large missiles set at squallx(squall,locust,hurricane), 
#large guns, and medium guns
analyze <- function(data) {

    #range with itu + bm + gi 
    #(+85 total for ballistic & energy, +0 for missiles) 
    gunsRangeMult <- 1.85
    
    mediumGuns <- read.csv("medium_guns.csv")
    largeGuns <- read.csv("large_guns.csv")
    mediumMissiles <- read.csv("medium_missiles.csv")
    largeMissiles <- read.csv("large_missiles.csv")
    
    dataFrame <- getDataFrame(data)
    dataFrameSorted <- getDataFrameSorted(data, dataFrame)
    
    #assemble tables
    shipTable <- getTKKTable(dataFrameSorted, Time ~ Ship, "Ship")             
    mediumGunTable <- getTTKTable(dataFrame, Time ~ mediumGuns, "Medium Gun")
    largeGunTable <- getTTKTable(dataFrame, Time ~ largeGuns, "Large Gun")
    mediumMissileTable <- getTTKTable(dataFrame, Time ~ mediumMissile, 
                                      "Medium Missile")
    largeMissileTable <- getTTKTable(dataFrame, Time ~ largeMissile, 
                                     "Large Missile")
    gunAndLargeMissilesTable <- getGunAndLargeMissileTable(dataFrameSorted)
    
    return(list(shipTable, "shipTable.txt",
                mediumGunTable, "mediumGunTable3.txt",
                largeGunTable, "largeGunTable3.txt",
                mediumMissileTable, "mediumMissileTable3.txt",
                largeMissileTable, "largeMissileTable3.txt",
                gunAndLargeMissilesTable, "gunAndLargeMissileTable.txt"))
}


main <- function() {
    tablesAndNames <- analyze(data)
    for (i in seq(1,tablesAndNames.length,2)) {
        writeTable(tablesAndNames[i], tablesAndNames[i+1]) 
    }
}

main()

[close]

CSVs

Code: medium_guns.csv

"Weapon Name","Flux Per Second","Range"
"Heavy Mortar",180,700
"Heavy Autocannon",214,800
"Arbalest",125,700
"Heavy Mauler",120,1000
"Hypervelocity Driver",175,1000
"Heavy Needler",200,700

Code: large_guns.csv

"Weapon Name","Flux Per Second","Range"
"Gauss",600,1200
"Mark IX",348,900
"Hellbore",250,900
"Storm Needler",350,700
"Mjolnir",667,900
"Hephaestus",480,900

Code: medium_missiles.csv

"Weapon Name","Flux Per Second","Range"
"Sabot",0,1400
"Harpoon",0,2500

Code: large_missiles.csv

"Weapon Name","Flux Per Second","Range"
"Squall",0,2500
"Locust",0,1400
"Hurricane",0,2500

[close]

[CapnHector]

Wow, you are fast and smart. Thanks for the effort and the contribution. It seems like this rewritten script does not actually write the necessary CSVs from text files, so that needs to be added.

Unfortunately some real life research work is becoming urgent. So despite being a little obsessed with Starsector just now I absolutely must set this aside for a couple of days at least and spend free time working on other things (I have somehow ended up trying to do a PhD while also doing two other degrees and working full time, the upside is there is never time to think about why you would ever want to do such a thing).

Anyway I think you mean you need the txt file (really tsv but I like Windows notepad) data containing the results from running the first script. So here it is. https://filetransfer.io/data-package/aWFSUaVb#link

I can't comment on every step of the code just now but here is a short version and I can go in more depth at another timepoint, also with some extra comments on already commented sections

Spoiler

# this is only needed for data visualization - can be safely removed
library(ggplot2)

#general functions
#engagementrange
range <- 1000
#fudge factor
errorsd <- 0.05
#the fudge factor should be a function of range (more error in position at greater range), but not a function of weapon firing angle, and be expressed in terms of pixels
error <- errorsd*range
#where does one shot hit within the weapon's arc of fire, in pixel difference from the target? get a random angle in degrees according to a uniform distribution, 
#then consider that the circumference is 2 pi * range pixels, so the hit coordinates in pixels are
shotangle <- function(acc) return(runif(1,-acc/2,acc/2)/360*2*pi*range)
#now add a random positional error to the coordinates of the hit
hitlocation <- function(acc){
  location <- shotangle(acc)
  location <- location + rnorm(1,0,error)
  return(location)
}
#so which box was hit?
cellhit <- function(angle){
  if(angle < anglerangevector[1]) return(1)
  if(angle > anglerangevector[ship[6]+1]) return(ship[6]+2)
  for (i in 1:length(anglerangevector)) {
    if ((angle > anglerangevector) & (angle <= anglerangevector[i+1])) return(i+1)
  }
}
# this function generates the shot distribution per 1 shot with 100000 samples
#note: it is likely that a sample size of 1000 or 10000 would be acceptable in practice, reducing distribution computation time meaningfully, I just wanted to go big
createdistribution <- function(acc){
  distributionvector <- vector(mode="double", length = ship[6]+2)
  for (i in 1:100000){
    wherehit <- cellhit(hitlocation(acc))
    distributionvector[wherehit] <- distributionvector[wherehit] +1
  }
  return(distributionvector/sum(distributionvector))
}
# this is the default distribution of damage to armor cells
b <- matrix(0,nrow=5,ncol=5)
b[1:5,2:4] <- 1/30
b[2:4,1:5] <- 1/30
b[2:4,2:4] <- 1/15
b[1,1] <- 0
b[1,5] <- 0
b[5,1] <- 0
b[5,5] <- 0
#this function generates a sum of matrices multiplied by the distribution

createhitmatrix <- function(acc){
  hitmatrix <- matrix(0,5,ship[6]+4)
  distributionvector <- createdistribution(acc)
  for (i in 1:ship[6]){
    hitmatrix[,i:(i+4)] <- hitmatrix[,i:(i+4)]+b*(distributionvector[i+1])
  }
  return(hitmatrix)
}

hitchance <- function(acc){
  hitchance <- 0
  distributionvector <- createdistribution(acc)
  hitchance <- (1-(distributionvector[1]+distributionvector[ship[6]]))
  return(hitchance)
}
#for weapons with damage changing over time we need a sequence of matrices
createhitmatrixsequence <- function(accvector){
  hitmatrixsequence <- list()
  for (i in 1:length(accvector)){
    hitmatrixsequence[] <- createhitmatrix(accvector)
  }
  return(hitmatrixsequence)
}
# we do not actually use the armordamage function in practice
armordamage <- function(damage, armor, startingarmor) damage*(max(0.15,damage/(damage+max(0.05*startingarmor,armor)))) 
# this atrociously named function contains a logical switch which probably degrades performance, but also ensures we never divide by 0
armordamageselectivereduction <- function(damage, armor,startingarmor) {
  useminarmor <- 0
  if(armor < 0.05*startingarmor / 15) useminarmor <- 1
  if(useminarmor == 0){
    if(armor == 0) {return (damage)}
    return(damage*(max(0.15,damage/(damage+armor))))
  }
  else{
    return(damage*(max(0.15,damage/(damage+0.05*startingarmor/15))))
  }
}


#ships - 
#ship, hullhp, shieldregen, shieldmax, startingarmor, widthinpixels, armorcells, name
glimmer <- c(1500, 250/0.6, 2500/0.6, 200, 78, 5, "glimmer")
brawlerlp <- c(2000, 500/0.8, 3000/0.8, 450,110,floor(110/15), "brawlerlp")
vanguard <- c(3000, 150, 2000, 600, 104, floor(104/15),"vanguard")
tempest <- c(1250, 225/0.6, 2500/0.6, 200,64,floor(64/15), "tempest")
medusa <- c(3000,400/0.6,6000/0.6,300,134,floor(134/15), "medusa")
hammerhead <- c(5000,250/0.8,4200/0.8,500,108,floor(108/16.4), "hammerhead")
enforcer <- c(4000,200,4000,900,136,floor(136/15), "enforcer")
dominator <- c(14000, 500, 10000, 1500, 180, 12, "dominator")
fulgent <- c(5000,300/0.6,5000/0.6,450, 160, floor(160/15), "fulgent")
brilliant <- c(8000,600/0.6,10000/0.6,900,160,floor(160/20),"brilliant")
radiant <- c(20000,1500/0.6,25000/0.6,1500,316,floor(316/30),"radiant")
onslaught <- c(20000,600,17000,1750,288,floor(288/30),"onslaught")
aurora <- c(8000,800/0.8,11000/0.8,800,128,floor(128/28), "aurora")
paragon <- c(18000,1250/0.6,25000/0.6,1500,330,floor(330/30),"paragon")
conquest <- c(12000,1200/1.4,20000/1.4,1200,190,floor(190/30),"conquest")
champion <- c(10000,550/0.8,10000/0.8,1250, 180,floor(180/24),"champion")

ships <- list(glimmer,brawlerlp,vanguard,tempest,medusa,hammerhead,enforcer,dominator,fulgent,brilliant,radiant,onslaught,aurora,paragon,conquest,champion)
#ships <- list(glimmer,brawlerlp)


#WEAPON ACCURACY
#missiles do not have spread
squallacc <- c(0)
locustacc <- c(0)
hurricaneacc <- c(0)
harpoonacc <- c(0)
sabotacc <- c(0)

#gauss has a spread of 0 and no increase per shot
gaussacc <- c(0)
#hephaestus has a spread of 0 and it increases by 2 per shot to a max of 10
hephaestusacc <- c(seq(0,10,2))
#mark ix has a spread of 0 and it increases by 2 per shot to a max of 15
markixacc <- c(seq(0,15,2),15)
#mjolnir has a spread of 0 and it increases by 1 per shot to a max of 5
mjolniracc <- c(seq(0,5,1))
#hellbore has a spread of 10
hellboreacc <- c(10)
#storm needler has a spread of 10
stormneedleracc <- c(10)

#squall fires 2 missiles / sec for 10 secs, then recharges for 10 secs
squalltics <- c(2,2,2,2,2,2,2,2,2,2,0,0,0,0,0,0,0,0,0,0)
#locust fires 10 missiles / sec for 4 secs, then recharges for 5 secs
locusttics <- c(10,10,10,10,0,0,0,0,0)
#hurricane fires 9 missiles every 15 seconds
hurricanetics <- c(9,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
#harpoon pod fires 4 missiles in 1 second, then recharges for 8 seconds (in reality 8.25)
harpoontics <- c(4,0,0,0,0,0,0,0,0)
#sabot pod fires 2*5 missiles, then recharges for 8 seconds (in reality 8.75 seconds and firing time is .5 seconds)
sabottics <- c(10,0,0,0,0,0,0,0,0)
#gauss fires 1 shot, then charges for 1 second
gausstics <- c(1,0)
#hephaestus fires 4 shots / sec
hephaestustics <- c(4)
#mark ix fires 4 shots every 3 seconds, taking .4 seconds for the burst, so 20 shots / 17 seconds in bursts of 4
markixtics <- c(4,0,0,4,0,0,4,0,0,4,0,0,0,4,0,0,0)
#mjolnir fires 4 shots per 3 seconds
mjolnirtics <- c(1,1,2,1,2,1,2,1,1)
#hellbore fires 1 shot per 4 seconds
hellboretics <- c(1,0,0,0)
#storm needler fires 10 shots per second
stormneedlertics <- c(10)

#arbalest fires 1 shot every 1.2 seconds, so 5 shots / 6 seconds
arbalesttics <- c(1,1,1,0,1,1)
arbalestacc <- c(0,5,10,15)
arbalest <- list(200, 2, arbalesttics, "Arbalest",arbalestacc)
#heavy mauler fires a burst of 3, taking .9 seconds to do so, then recharges for 4.5 seconds. So total of 3 shots / 5.4 seconds, so total of 15 shots / 27 seconds in bursts of 3.
heavymaulertics <- c(3,0,0,0,0,0,3,0,0,0,0,3,0,0,0,0,3,0,0,0,0,3,0,0,0,0,0)
heavymauleracc <- c(seq(0,5,1))
heavymauler <- list(200,0.5,heavymaulertics, "Heavy Mauler", heavymauleracc)
#heavy autocannon takes .6 seconds to fire a burst, then recharges for 1 second. So 3 shots / 1.6 secs, so 15 shots / 8 secs in bursts
hactics <- c(3,0,3,3,0,3,3,0)
hacacc <- c(seq(0,18,3))
hac <- list(100, 2, hactics, "Heavy Autocannon", hacacc)
hvdtics <- c(1,0,0)
hvdacc <- c(0)
hvd <- list(275,2,hvdtics,"Hypervelocity Driver", hvdacc)
#heavy needler fires a burst of 30 needles with 0.05 sec intervals (1.5 sec) and then recharges for 4.55 sec (round to 4.5 sec)
needlertics <- c(20,10,0,0,0,0)
needleracc <- c(seq(1,10,0.5))
needler <- list(50,2,needlertics,"Heavy Needler",needleracc)
#heavy mortar fires a burst of 2 rounds, taking .6 seconds to do so, then recharges for .7 seconds. So a total of 2 rounds every 1.3 seconds, so 20 rounds / 13 seconds in bursts of 2.
mortartics <- c(2,2,0,2,2,2,2,0,2,2,2,0,2)
mortaracc <- c(seq(0,20,5))
mortar <- list(110,0.5,mortartics, "Heavy Mortar", mortaracc)

#damage per shot, damage type (2=kinetic, 0.5=he, 0.25=frag, 1=energy), tics, weapon name, weapon accuracy over time, hit chance
squall <- list(250, 2, squalltics, "Squall", squallacc)
locust <- list(200, 0.25, locusttics, "Locust", locustacc)
hurricane <- list(500, 0.5, hurricanetics, "Hurricane", hurricaneacc)
harpoon <- list(750, 0.5, harpoontics, "Harpoon", harpoonacc)
sabot <- list(200, 2, sabottics, "Sabot", sabotacc)
gauss <- list(700, 2, gausstics, "Gauss", gaussacc)
hephaestus <- list(120, 0.5, hephaestustics, "Hephaestus", hephaestusacc)
markix <- list(200, 2, markixtics, "Mark IX", markixacc)
mjolnir <- list(400, 1, mjolnirtics, "Mjolnir", mjolniracc)
hellbore <- list(750, 0.5, hellboretics, "Hellbore", hellboreacc)
stormneedler <- list(50, 2, stormneedlertics, "Storm Needler", stormneedleracc)
dummy <- list(0,0,c(0),"",c(0),c(0))

#which weapons are we studying?

weapon1choices <- list(squall, locust, hurricane)
weapon2choices <- list(squall, locust, hurricane)
weapon3choices <- list(harpoon, sabot)
weapon4choices <- list(harpoon, sabot)
weapon5choices <- list(gauss, markix, mjolnir, hellbore, hephaestus, stormneedler)
weapon6choices <- list(gauss, markix, mjolnir, hellbore, hephaestus, stormneedler)
weapon7choices <- list(arbalest, hac, hvd, heavymauler, needler, mortar)
weapon8choices <- list(arbalest, hac, hvd, heavymauler, needler, mortar)

#how many unique weapon loadouts are there?

#get names of weapons from a choices list x
getweaponnames <- function(x){
  vector <- vector(mode="character")
  for (i in 1:length(x)){
  vector <- cbind(vector, x[][[4]])
  }
  return(vector)
}
#convert the names back to numbers when we are done based on a weapon choices list y
convertweaponnames <- function(x, y){
  vector <- vector(mode="integer")
  for (j in 1:length(x)) {
    for (i in 1:length(y)){
      if(x[j] == y[][[4]]) vector <- cbind(vector, i)
    }
  }
  return(vector)
}

#this section of code generates a table of all unique loadouts that we can create using the weapon choices available
generatepermutations <- 1
if (generatepermutations == 1){
#enumerate weapon choices as integers

perm1 <- seq(1,length(weapon1choices),1)
perm2 <- seq(1,length(weapon2choices),1)
perm3 <- seq(1,length(weapon3choices),1)
perm4 <- seq(1,length(weapon4choices),1)
perm5 <- seq(1,length(weapon5choices),1)
perm6 <- seq(1,length(weapon6choices),1)
perm7 <- seq(1,length(weapon7choices),1)
perm8 <- seq(1,length(weapon8choices),1)

#create a matrix of all combinations
perm1x2 <- expand.grid(perm1,perm2)
#sort, then only keep unique rows
perm1x2 <- unique(t(apply(perm1x2, 1, sort)))

perm3x4 <- expand.grid(perm3,perm4)
perm3x4 <- unique(t(apply(perm3x4, 1, sort)))

perm5x6 <- expand.grid(perm5,perm6)
perm5x6 <- unique(t(apply(perm5x6, 1, sort)))

perm7x8 <- expand.grid(perm7,perm8)
perm7x8 <- unique(t(apply(perm7x8, 1, sort)))

#now that we have all unique combinations of all two weapons, create a matrix containing all combinations of these unique combinations
allperms <- matrix(0,0,(length(perm1x2[1,])+length(perm3x4[1,])+length(perm5x6[1,])+length(perm7x8[1,])))
for(i in 1:length(perm1x2[,1])) for(j in 1:length(perm3x4[,1])) for(k in 1:length(perm5x6[,1])) for(l in 1:length(perm7x8[,1])) allperms <- rbind(allperms, c(perm1x2[i,],perm3x4[j,],perm5x6[k,],perm7x8[l,])
)
#this is just for testing, can remove
allperms
#we save this so we don't have to compute it again
saveRDS(allperms, file="allperms.RData")

} else {
  allperms <- readRDS("lookuptable.RData")
}

#now compute a main lookuptable to save on computing time
#the lookuptable should be a list of lists, so that 
#lookuptable[[ship]][[weapon]][[1]] returns hit chance vector and 
#lookuptable[[ship]][[weapon]][[2]] returns hit probability matrix
#time for some black R magic

#note: the lookuptable will be formulated such that there is a running index of weapons rather than sub-lists, so all weapons will be indexed consecutively so we have lookuptable [[1]][[1]] = [[ship1]][[weaponchoices1_choice1]], etc. So that is what the below section does.

#read or generate lookuptable
generatelookuptable <- 0
if(generatelookuptable == 1){

lookuptable <- list()

for (f in 1:length(ships)){
  lookuptable[[f]] <- list()
  ship <- ships[[f]]
  ship <- as.double(ship[1:6])
  #how much is the visual arc of the ship in rad?
  shipangle <- ship[5]/(2* pi *range)
  
  #how much is the visual arc of a single cell of armor in rad?
  cellangle <- shipangle/ship[6]
  
  #now assume the weapon is targeting the center of the ship's visual arc and that the ship is in the center of the weapon's firing arc
  #which cell will the shot hit, or will it miss?
  #call the cells (MISS, cell1, cell2, ... ,celli, MISS) and get a vector giving the (maximum for negative / minimum for positive) angles for hitting each
  anglerangevector <- vector(mode="double", length = ship[6]+1)
  anglerangevector[1] <- -shipangle/2
  for (i in 1:(length(anglerangevector)-1)) anglerangevector[i+1] <- anglerangevector+cellangle
  #now convert it to pixels
  anglerangevector <- anglerangevector*2*pi*range
  
weaponindexmax <- length(weapon1choices)+length(weapon2choices)+length(weapon3choices)+length(weapon4choices)+length(weapon5choices)+length(weapon6choices)+length(weapon7choices)+length(weapon8choices)
  
for (x in 1:weaponindexmax) {
  print(x)
  if(x <= length(weapon1choices)){
  weapon1<-weapon1choices[
]

  if(weapon1[4] != ""){
    hitchancevector <- vector(mode = "double", length = length(weapon1[[5]]))
    for (i in 1:length(weapon1[[5]])){
      hitchancevector <- hitchance(weapon1[[5]])
    }
    lookuptable[[f]][
] <- list()

    lookuptable[[f]][
][[1]] <- hitchancevector

    lookuptable[[f]][
][[2]] <- createhitmatrixsequence(weapon1[[5]])

  }
  }
  if((x > length(weapon1choices)) &  (x <= length(weapon1choices) + length(weapon2choices))){
    weapon2<-weapon2choices[[x-length(weapon1choices)]]
    if(weapon2[4] != ""){
      hitchancevector <- vector(mode = "double", length = length(weapon2[[5]]))
      for (i in 1:length(weapon2[[5]])){
        hitchancevector <- hitchance(weapon2[[5]])
      }
      lookuptable[[f]][
] <- list()

      lookuptable[[f]][
][[1]] <- hitchancevector

      lookuptable[[f]][
][[2]] <- createhitmatrixsequence(weapon2[[5]])

    }
  }
  
  if((x > length(weapon1choices) + length(weapon2choices)) &  (x <= length(weapon2choices) + length(weapon1choices) + length(weapon3choices))){
    weapon3<-weapon3choices[[x-length(weapon2choices)-length(weapon1choices)]]
    if(weapon3[4] != ""){
      hitchancevector <- vector(mode = "double", length = length(weapon3[[5]]))
      for (i in 1:length(weapon3[[5]])){
        hitchancevector <- hitchance(weapon3[[5]])
      }
      lookuptable[[f]][
] <- list()

      lookuptable[[f]][
][[1]] <- hitchancevector

      lookuptable[[f]][
][[2]] <- createhitmatrixsequence(weapon3[[5]])

    }
  }  
  
  if((x > length(weapon2choices) + length(weapon1choices) + length(weapon3choices)) &  (x <= length(weapon3choices) + length(weapon2choices) + length(weapon1choices) + length(weapon4choices))){
    weapon4<-weapon4choices[[x-length(weapon3choices)-length(weapon2choices)-length(weapon1choices)]]
    if(weapon4[4] != ""){
      hitchancevector <- vector(mode = "double", length = length(weapon4[[5]]))
      for (i in 1:length(weapon4[[5]])){
        hitchancevector <- hitchance(weapon4[[5]])
      }
      lookuptable[[f]][
] <- list()

      lookuptable[[f]][
][[1]] <- hitchancevector

      lookuptable[[f]][
][[2]] <- createhitmatrixsequence(weapon4[[5]])

    }
  }
  if((x > length(weapon3choices) + length(weapon2choices) + length(weapon1choices) + length(weapon4choices)) &  (x <= length(weapon4choices) + length(weapon3choices) + length(weapon2choices) + length(weapon1choices) + length(weapon5choices))){
    weapon5<-weapon5choices[[x-length(weapon4choices)-length(weapon3choices)-length(weapon2choices)-length(weapon1choices)]]
    if(weapon5[4] != ""){
      hitchancevector <- vector(mode = "double", length = length(weapon5[[5]]))
      for (i in 1:length(weapon5[[5]])){
        hitchancevector <- hitchance(weapon5[[5]])
      }
      lookuptable[[f]][
] <- list()

      lookuptable[[f]][
][[1]] <- hitchancevector

      lookuptable[[f]][
][[2]] <- createhitmatrixsequence(weapon5[[5]])

    }
  }
  if((x > length(weapon4choices) + length(weapon3choices) + length(weapon2choices) + length(weapon1choices) + length(weapon5choices)) &  (x <= length(weapon5choices) + length(weapon4choices) + length(weapon3choices) + length(weapon2choices) + length(weapon1choices) + length(weapon6choices))){
    weapon6<-weapon6choices[[x-length(weapon5choices)-length(weapon4choices)-length(weapon3choices)-length(weapon2choices)-length(weapon1choices)]]
    if(weapon6[4] != ""){
      hitchancevector <- vector(mode = "double", length = length(weapon6[[5]]))
      for (i in 1:length(weapon6[[5]])){
        hitchancevector <- hitchance(weapon6[[5]])
      }
      lookuptable[[f]][
] <- list()

      lookuptable[[f]][
][[1]] <- hitchancevector

      lookuptable[[f]][
][[2]] <- createhitmatrixsequence(weapon6[[5]])

    }
  }
  if((x > length(weapon5choices) + length(weapon4choices) + length(weapon3choices) + length(weapon2choices) + length(weapon1choices) + length(weapon6choices)) &  (x <= length(weapon6choices) + length(weapon5choices) + length(weapon4choices) + length(weapon3choices) + length(weapon2choices) + length(weapon1choices) + length(weapon7choices))){
    weapon7<-weapon7choices[[x-length(weapon6choices)-length(weapon5choices)-length(weapon4choices)-length(weapon3choices)-length(weapon2choices)-length(weapon1choices)]]
    if(weapon7[4] != ""){
      hitchancevector <- vector(mode = "double", length = length(weapon7[[5]]))
      for (i in 1:length(weapon7[[5]])){
        hitchancevector <- hitchance(weapon7[[5]])
      }
      lookuptable[[f]][
] <- list()

      lookuptable[[f]][
][[1]] <- hitchancevector

      lookuptable[[f]][
][[2]] <- createhitmatrixsequence(weapon7[[5]])

    }
  }
  if((x > length(weapon6choices) + length(weapon5choices) + length(weapon4choices) + length(weapon3choices) + length(weapon2choices) + length(weapon1choices) + length(weapon7choices)) &  (x <= length(weapon7choices) + length(weapon6choices) + length(weapon5choices) + length(weapon4choices) + length(weapon3choices) + length(weapon2choices) + length(weapon1choices) + length(weapon8choices))){
    weapon8<-weapon8choices[[x-length(weapon7choices)-length(weapon6choices)-length(weapon5choices)-length(weapon4choices)-length(weapon3choices)-length(weapon2choices)-length(weapon1choices)]]
    if(weapon8[4] != ""){
      hitchancevector <- vector(mode = "double", length = length(weapon8[[5]]))
      for (i in 1:length(weapon8[[5]])){
        hitchancevector <- hitchance(weapon8[[5]])
      }
      lookuptable[[f]][
] <- list()

      lookuptable[[f]][
][[1]] <- hitchancevector

      lookuptable[[f]][
][[2]] <- createhitmatrixsequence(weapon8[[5]])

    }
  }
  
}
}
# save this so we don't have to re-compute it
saveRDS(lookuptable, file="lookuptable.RData")
} else {
lookuptable <- readRDS("lookuptable.RData")
}

lookup <- function(ship, weapon, var) return(lookuptable[[ship]][[weapon]][[var]])
                    
#this generates a heatmap for visualization, remove
geom_tile(lookup(8,15,2)[[2]])

#go through all ships
for (f in 1:length(ships)){

ship <- ships[[f]]
#format ship data types appropriately
shipname <- ship[[7]]
ship <- as.double(ship[1:6])







timeseriesarray <- data.frame(matrix(ncol = 7,nrow=0))



timetokill=0


shieldblock <- 0

#calculate shield damage
#timepoint %% length ... etc. section returns how many shots that weapon will fire at that point of it's cycle (ie. vector index = timepoint modulo vector index maximum)

shielddamageattimepoint <- function(weapon, timepoint){
  nohits <- weapon[[3]][(timepoint %% (length(weapon[[3]])))+1]
  if (nohits == 0) {return(0)} else {
    return(weapon[[1]]*nohits)
  }
}

#hitmatrix[hitx,hity] returns the damage allocation to that cell

damageattimepoint <- function(weapon, timepoint, armor, startingarmor, hitx, hity, shots){
  # vectors in R are indexed starting from 1
  hitmatrix <- weapon[[7]][[min(shots,length(weapon[[7]]))]]
  nohits <- weapon[[3]][(timepoint %% (length(weapon[[3]])))+1]
  if (nohits == 0) {return(0)} else {
    damagesum <- 0
    for (i in 1:nohits) {
      #      damagesum <- damagesum + as.double(weapon[[1]]*hitmatrix[hitx,hity])
      damagesum <- damagesum + armordamageselectivereduction(as.double(weapon[[1]]*hitmatrix[hitx,hity]),armor,startingarmor)
      shots <- shots + 1
      hitmatrix <- weapon[[7]][[min(shots,length(weapon[[7]]))]]
    }
    return(damagesum)
  }
}



timeseries <- function(timepoint, shieldhp, armorhp, hullhp, shieldregen, shieldmax, startingarmor,armormatrix){
  weaponacc <- 0
  #are we using shield to block?
  shieldblock <- 0
  hulldamage <- 0
  if(hullhp > 0){} else {shieldhp <- 0}
  
  #1. shields. if shieldhp is sufficient, use shield to block
  if (weapon1[[4]] !=""){
    weapon1mult <- unlist(weapon1[2])
    if (shieldhp > shielddamageattimepoint(weapon1, timepoint)*weapon1mult){
      shieldhp <- shieldhp - shielddamageattimepoint(weapon1, timepoint)*weapon1mult*weapon1[[6]][min(weapon1shots,length(weapon1[[6]]))]
      shieldhp <- max(shieldhp, 0)
      if(shielddamageattimepoint(weapon1,timepoint) > 0) {shieldblock <- 1}
    } else {
      #if you did not use shield to block, regenerate flux <- applied later
      #2. armor and hull
      if(unlist(weapon1[2])==0.25){weapon1mult = 0.25} else {weapon1mult= 1 / unlist(weapon1[2])}
      #2.1. damage armor and hull
      hulldamage <- 0    
      
      #damageattimepoint <- function(weapon, timepoint, armor, startingarmor, hitx, hity, shots){
      
      for (j in 1:length(armormatrix[1,])){
        for (i in 1:length(armormatrix[,1])){

#this monster of a line does the following:
#hull damage is maximum of: 0, or (weapon damage to armor cell, with multiplier - armor cell hp)/weapon multiplier
          hulldamage <- hulldamage+max(0,weapon1mult*damageattimepoint(weapon1, timepoint,armormatrix[i,j],startingarmor,i,j,weapon1shots)-armormatrix[i,j])/weapon1mult
#new armor hp at cell [i,j] is maximum of: 0, or (armor cell health - weapon damage to that section of armor at that timepoint, multiplied by overall hit probability to hit ship --- note: the damage distribution matrix is calculated such that it is NOT inclusive of hit chance, but hit chance is added separately here
          armormatrix[i,j] <- max(0,armormatrix[i,j]-weapon1mult*damageattimepoint(weapon1, timepoint,armormatrix[i,j],startingarmor,i,j,weapon1shots)*weapon1[[6]][min(weapon1shots,length(weapon1[[6]]))])
        }}
      
      # now reduce hullhp by hull damage times probability to hit ship
      hullhp <- hullhp - hulldamage*weapon1[[6]][min(weapon1shots,length(weapon1[[6]]))]
      hullhp <- max(hullhp, 0)
      
      
    }
  }
  weapon1shots <- weapon1shots + weapon1[[3]][(timepoint %% (length(weapon1[[3]])+1))]
  
  
  #repeat for other weapons
  #1. shields. if shieldhp is sufficient, use shield to block
  if (weapon2[[4]] !=""){
    weapon2mult <- unlist(weapon2[2])
    if (shieldhp > shielddamageattimepoint(weapon2, timepoint)*weapon2mult){
      shieldhp <- shieldhp - shielddamageattimepoint(weapon2, timepoint)*weapon2mult*weapon2[[6]][min(weapon2shots,length(weapon2[[6]]))]
      shieldhp <- max(shieldhp, 0)
      if(shielddamageattimepoint(weapon2,timepoint) > 0) {shieldblock <- 1}
    } else {
      #if you did not use shield to block, regenerate flux
      #2. armor and hull
      if(unlist(weapon2[2])==0.25){weapon2mult = 0.25} else {weapon2mult= 1 / unlist(weapon2[2])}
      #2.1. damage armor and hull
      hulldamage <- 0    
      
      #damageattimepoint <- function(weapon, timepoint, armor, startingarmor, hitx, hity, shots){
      
      for (j in 1:length(armormatrix[1,])){
        for (i in 1:length(armormatrix[,1])){
          hulldamage <- hulldamage+max(0,weapon2mult*damageattimepoint(weapon2, timepoint,armormatrix[i,j],startingarmor,i,j,weapon2shots)-armormatrix[i,j])/weapon2mult
          armormatrix[i,j] <- max(0,armormatrix[i,j]-weapon2mult*damageattimepoint(weapon2, timepoint,armormatrix[i,j],startingarmor,i,j,weapon2shots)*weapon2[[6]][min(weapon2shots,length(weapon2[[6]]))])
        }}
      
      
      hullhp <- hullhp - hulldamage*weapon2[[6]][min(weapon2shots,length(weapon2[[6]]))]
      hullhp <- max(hullhp, 0)
      
      
    }
  }
  weapon2shots <- weapon2shots + weapon2[[3]][(timepoint %% (length(weapon2[[3]])+1))]
  #1. shields. if shieldhp is sufficient, use shield to block
  if (weapon3[[4]] !=""){
    weapon3mult <- unlist(weapon3[2])
    if (shieldhp > shielddamageattimepoint(weapon3, timepoint)*weapon3mult){
      shieldhp <- shieldhp - shielddamageattimepoint(weapon3, timepoint)*weapon3mult*weapon3[[6]][min(weapon3shots,length(weapon3[[6]]))]
      shieldhp <- max(shieldhp, 0)
      if(shielddamageattimepoint(weapon3,timepoint) > 0) {shieldblock <- 1}
    } else {
      #if you did not use shield to block, regenerate flux
      #2. armor and hull
      if(unlist(weapon3[2])==0.25){weapon3mult = 0.25} else {weapon3mult= 1 / unlist(weapon3[2])}
      #2.1. damage armor and hull
      hulldamage <- 0    
      
      #damageattimepoint <- function(weapon, timepoint, armor, startingarmor, hitx, hity, shots){
      
      for (j in 1:length(armormatrix[1,])){
        for (i in 1:length(armormatrix[,1])){
          hulldamage <- hulldamage+max(0,weapon3mult*damageattimepoint(weapon3, timepoint,armormatrix[i,j],startingarmor,i,j,weapon3shots)-armormatrix[i,j])/weapon3mult
          armormatrix[i,j] <- max(0,armormatrix[i,j]-weapon3mult*damageattimepoint(weapon3, timepoint,armormatrix[i,j],startingarmor,i,j,weapon3shots)*weapon3[[6]][min(weapon3shots,length(weapon3[[6]]))])
        }}
      
      
      hullhp <- hullhp - hulldamage*weapon3[[6]][min(weapon3shots,length(weapon3[[6]]))]
      hullhp <- max(hullhp, 0)
      
      
    }
  }
  weapon3shots <- weapon3shots + weapon3[[3]][(timepoint %% (length(weapon3[[3]])+1))]
  #1. shields. if shieldhp is sufficient, use shield to block
  if (weapon4[[4]] !=""){
    weapon4mult <- unlist(weapon4[2])
    if (shieldhp > shielddamageattimepoint(weapon4, timepoint)*weapon4mult){
      shieldhp <- shieldhp - shielddamageattimepoint(weapon4, timepoint)*weapon4mult*weapon4[[6]][min(weapon4shots,length(weapon4[[6]]))]
      shieldhp <- max(shieldhp, 0)
      if(shielddamageattimepoint(weapon4,timepoint) > 0) {shieldblock <- 1}
    } else {
      #if you did not use shield to block, regenerate flux
      #2. armor and hull
      if(unlist(weapon4[2])==0.25){weapon4mult = 0.25} else {weapon4mult= 1 / unlist(weapon4[2])}
      #2.1. damage armor and hull
      hulldamage <- 0    
      
      #damageattimepoint <- function(weapon, timepoint, armor, startingarmor, hitx, hity, shots){
      
      for (j in 1:length(armormatrix[1,])){
        for (i in 1:length(armormatrix[,1])){
          hulldamage <- hulldamage+max(0,weapon4mult*damageattimepoint(weapon4, timepoint,armormatrix[i,j],startingarmor,i,j,weapon4shots)-armormatrix[i,j])/weapon4mult
          armormatrix[i,j] <- max(0,armormatrix[i,j]-weapon4mult*damageattimepoint(weapon4, timepoint,armormatrix[i,j],startingarmor,i,j,weapon4shots)*weapon4[[6]][min(weapon4shots,length(weapon4[[6]]))])
        }}
      
      
      hullhp <- hullhp - hulldamage*weapon4[[6]][min(weapon4shots,length(weapon4[[6]]))]
      hullhp <- max(hullhp, 0)
      
      
    }
  }
  weapon4shots <- weapon4shots + weapon4[[3]][(timepoint %% (length(weapon4[[3]])+1))]
  #1. shields. if shieldhp is sufficient, use shield to block
  if (weapon5[[4]] !=""){
    weapon5mult <- unlist(weapon5[2])
    if (shieldhp > shielddamageattimepoint(weapon5, timepoint)*weapon5mult){
      shieldhp <- shieldhp - shielddamageattimepoint(weapon5, timepoint)*weapon5mult*weapon5[[6]][min(weapon5shots,length(weapon5[[6]]))]
      shieldhp <- max(shieldhp, 0)
      if(shielddamageattimepoint(weapon5,timepoint) > 0) {shieldblock <- 1}
    } else {
      #if you did not use shield to block, regenerate flux
      #2. armor and hull
      if(unlist(weapon5[2])==0.25){weapon5mult = 0.25} else {weapon5mult= 1 / unlist(weapon5[2])}
      #2.1. damage armor and hull
      hulldamage <- 0    
      
      #damageattimepoint <- function(weapon, timepoint, armor, startingarmor, hitx, hity, shots){
      
      for (j in 1:length(armormatrix[1,])){
        for (i in 1:length(armormatrix[,1])){
          hulldamage <- hulldamage+max(0,weapon5mult*damageattimepoint(weapon5, timepoint,armormatrix[i,j],startingarmor,i,j,weapon5shots)-armormatrix[i,j])/weapon5mult
          armormatrix[i,j] <- max(0,armormatrix[i,j]-weapon5mult*damageattimepoint(weapon5, timepoint,armormatrix[i,j],startingarmor,i,j,weapon5shots)*weapon5[[6]][min(weapon5shots,length(weapon5[[6]]))])
        }}
      
      
      hullhp <- hullhp - hulldamage*weapon5[[6]][min(weapon5shots,length(weapon5[[6]]))]
      hullhp <- max(hullhp, 0)
      
      
    }
  }
  weapon5shots <- weapon5shots + weapon5[[3]][(timepoint %% (length(weapon5[[3]])+1))]
  #1. shields. if shieldhp is sufficient, use shield to block
  if (weapon6[[4]] !=""){
    weapon6mult <- unlist(weapon6[2])
    if (shieldhp > shielddamageattimepoint(weapon6, timepoint)*weapon6mult){
      shieldhp <- shieldhp - shielddamageattimepoint(weapon6, timepoint)*weapon6mult*weapon6[[6]][min(weapon6shots,length(weapon6[[6]]))]
      shieldhp <- max(shieldhp, 0)
      if(shielddamageattimepoint(weapon6,timepoint) > 0) {shieldblock <- 1}
    } else {
      #if you did not use shield to block, regenerate flux
      #2. armor and hull
      if(unlist(weapon6[2])==0.25){weapon6mult = 0.25} else {weapon6mult= 1 / unlist(weapon6[2])}
      #2.1. damage armor and hull
      hulldamage <- 0    
      
      #damageattimepoint <- function(weapon, timepoint, armor, startingarmor, hitx, hity, shots){
      
      for (j in 1:length(armormatrix[1,])){
        for (i in 1:length(armormatrix[,1])){
          hulldamage <- hulldamage+max(0,weapon6mult*damageattimepoint(weapon6, timepoint,armormatrix[i,j],startingarmor,i,j,weapon6shots)-armormatrix[i,j])/weapon6mult
          armormatrix[i,j] <- max(0,armormatrix[i,j]-weapon6mult*damageattimepoint(weapon6, timepoint,armormatrix[i,j],startingarmor,i,j,weapon6shots)*weapon6[[6]][min(weapon6shots,length(weapon6[[6]]))])
        }}
      
      
      hullhp <- hullhp - hulldamage*weapon6[[6]][min(weapon6shots,length(weapon6[[6]]))]
      hullhp <- max(hullhp, 0)
      
      
    }
  }
  weapon6shots <- weapon6shots + weapon6[[3]][(timepoint %% (length(weapon6[[3]])+1))]
  
  #1. shields. if shieldhp is sufficient, use shield to block
  if (weapon7[[4]] !=""){
    weapon7mult <- unlist(weapon7[2])
    if (shieldhp > shielddamageattimepoint(weapon7, timepoint)*weapon7mult){
      shieldhp <- shieldhp - shielddamageattimepoint(weapon7, timepoint)*weapon7mult*weapon7[[6]][min(weapon7shots,length(weapon7[[6]]))]
      shieldhp <- max(shieldhp, 0)
      if(shielddamageattimepoint(weapon7,timepoint) > 0) {shieldblock <- 1}
    } else {
      #if you did not use shield to block, regenerate flux
      #2. armor and hull
      if(unlist(weapon7[2])==0.25){weapon7mult = 0.25} else {weapon7mult= 1 / unlist(weapon7[2])}
      #2.1. damage armor and hull
      hulldamage <- 0    
      
      #damageattimepoint <- function(weapon, timepoint, armor, startingarmor, hitx, hity, shots){
      
      for (j in 1:length(armormatrix[1,])){
        for (i in 1:length(armormatrix[,1])){
          hulldamage <- hulldamage+max(0,weapon7mult*damageattimepoint(weapon7, timepoint,armormatrix[i,j],startingarmor,i,j,weapon7shots)-armormatrix[i,j])/weapon7mult
          armormatrix[i,j] <- max(0,armormatrix[i,j]-weapon7mult*damageattimepoint(weapon7, timepoint,armormatrix[i,j],startingarmor,i,j,weapon7shots)*weapon7[[6]][min(weapon7shots,length(weapon7[[6]]))])
        }}
      
      
      hullhp <- hullhp - hulldamage*weapon7[[6]][min(weapon7shots,length(weapon7[[6]]))]
      hullhp <- max(hullhp, 0)
      
      
    }
  }
  weapon7shots <- weapon7shots + weapon7[[3]][(timepoint %% (length(weapon7[[3]])+1))]
  
  #1. shields. if shieldhp is sufficient, use shield to block
  if (weapon8[[4]] !=""){
    weapon8mult <- unlist(weapon8[2])
    if (shieldhp > shielddamageattimepoint(weapon8, timepoint)*weapon8mult){
      shieldhp <- shieldhp - shielddamageattimepoint(weapon8, timepoint)*weapon8mult*weapon8[[6]][min(weapon8shots,length(weapon8[[6]]))]
      shieldhp <- max(shieldhp, 0)
      if(shielddamageattimepoint(weapon8,timepoint) > 0) {shieldblock <- 1}
    } else {
      #if you did not use shield to block, regenerate flux
      #2. armor and hull
      if(unlist(weapon8[2])==0.25){weapon8mult = 0.25} else {weapon8mult= 1 / unlist(weapon8[2])}
      #2.1. damage armor and hull
      hulldamage <- 0    
      
      #damageattimepoint <- function(weapon, timepoint, armor, startingarmor, hitx, hity, shots){
      
      for (j in 1:length(armormatrix[1,])){
        for (i in 1:length(armormatrix[,1])){
          hulldamage <- hulldamage+max(0,weapon8mult*damageattimepoint(weapon8, timepoint,armormatrix[i,j],startingarmor,i,j,weapon8shots)-armormatrix[i,j])/weapon8mult
          armormatrix[i,j] <- max(0,armormatrix[i,j]-weapon8mult*damageattimepoint(weapon8, timepoint,armormatrix[i,j],startingarmor,i,j,weapon8shots)*weapon8[[6]][min(weapon8shots,length(weapon8[[6]]))])
        }}
      
      
      hullhp <- hullhp - hulldamage*weapon8[[6]][min(weapon8shots,length(weapon8[[6]]))]
      hullhp <- max(hullhp, 0)
      
      
    }
  }
  weapon8shots <- weapon8shots + weapon8[[3]][(timepoint %% (length(weapon8[[3]])+1))]
  
  
  
  hullhp <- hullhp - hulldamage
  hullhp <- max(hullhp, 0)
  
  if(hullhp==0) armorhp <- 0
  
  
  
  
  armorhp <- sum(armormatrix)*15/((ship[[6]]+4)*5)
  if(hullhp==0) armorhp <- 0
  
  
  if (shieldblock==0){shieldhp <- min(shieldmax,shieldhp+shieldregen)}
  return(list(timepoint, shieldhp, armorhp, hullhp, shieldregen, shieldmax, startingarmor,armormatrix))
}

totaltime = 500



armorhp <- ship[4]
shieldhp <- ship[3]
hullhp <- ship[1]
shieldregen <- ship[2]
shieldmax <- ship[3]
armorhp <- ship[4]
startingarmor <- ship[4]
weapon1shots <- 1
weapon2shots <- 1
weapon3shots <- 1
weapon4shots <- 1
weapon5shots <- 1
weapon6shots <- 1
weapon7shots <- 1
weapon8shots <- 1

armormatrix <- matrix(ship[4]/15,5,ship[6]+4)
startingarmormatrixsum <- sum(armormatrix)

#now what we do here is we go through all the permutations using the running index, which is i+j+k+l+m+n+o+p for weapons 8
for (z in 1:length(allperms[,1])) {
  i <- allperms[z,1]
  j <- allperms[z,2]
  k <- allperms[z,3]
  l <- allperms[z,4]
  m <- allperms[z,5]
  n <- allperms[z,6]
  o <- allperms[z,7]
  p <- allperms[z,8]
  
#for (i in 1:length(weapon1choices)) {
  weapon1<-weapon1choices[]
#  for (j in 1:length(weapon2choices)) {
    weapon2<-weapon2choices[[j]]
#    for (k in 1:length(weapon3choices)) {
      weapon3<-weapon3choices[[k]]
#      for (l in 1:length(weapon4choices)) {
        weapon4<-weapon4choices[[l]]
#        for (m in 1:length(weapon5choices)) {
          weapon5<-weapon5choices[[m]]
#          for (n in 1:length(weapon6choices)) {
            weapon6<-weapon6choices[[n]]
#            for (o in 1:length(weapon7choices)) {
              weapon7<-weapon7choices[
]

#              for (p in 1:length(weapon8choices)) {
                weapon8<-weapon8choices[[p]]
                #lookup <- function(ship, weapon, var) return(lookuptable[[ship]][[weapon]][[var]])
                if(weapon1[4] != ""){
                  weapon1[[6]] <- lookup(f,i,1)
                  weapon1[[7]] <- lookup(f,i,2)
                }
                
                if(weapon2[4] != ""){
                  weapon2[[6]] <- lookup(f,i+j,1)
                  weapon2[[7]] <- lookup(f,i+j,2)
                }
                
                if(weapon3[4] != ""){
                  weapon3[[6]] <- lookup(f,i+j+k,1)
                  weapon3[[7]] <- lookup(f,i+j+k,2)
                }
                
                if(weapon4[4] != ""){
                  weapon4[[6]] <- lookup(f,i+j+k+l,1)
                  weapon4[[7]] <- lookup(f,i+j+k+l,2)
                }
                
                if(weapon5[4] != ""){
                  weapon5[[6]] <- lookup(f,i+j+k+l+m,1)
                  weapon5[[7]] <- lookup(f,i+j+k+l+m,2)
                }
                
                if(weapon6[4] != ""){
                  weapon6[[6]] <- lookup(f,i+j+k+l+m+n,1)
                  weapon6[[7]] <- lookup(f,i+j+k+l+m+n,2)
                }
                if(weapon7[4] != ""){
                  weapon7[[6]] <- lookup(f,i+j+k+l+m+n+o,1)
                  weapon7[[7]] <- lookup(f,i+j+k+l+m+n+o,2)
                }
                if(weapon8[4] != ""){
                  weapon8[[6]] <- lookup(f,i+j+k+l+m+n+o+p,1)
                  weapon8[[7]] <- lookup(f,i+j+k+l+m+n+o+p,2)
                }
                
#time series - run time series at point t, save it to state, update values according to state, re-run time series, break if ship dies
for (t in 1:totaltime){
  state <- timeseries(t,shieldhp,armorhp,hullhp,shieldregen,shieldmax,startingarmor,armormatrix)
  shieldhp <- state[[2]]
  armorhp <- state[[3]]
  hullhp <- state[[4]]
  flux <- shieldmax - shieldhp
  armormatrix <- state[[8]]
  if(hullhp == 0){flux <- 0
  if (timetokill == 0){timetokill <- t
  break}
  }
  
}
if (timetokill ==0){timetokill <- NA}

tobind <- c(timetokill,unlist(weapon1[4]),unlist(weapon2[4]),unlist(weapon3[4]),unlist(weapon4[4]),unlist(weapon5[4]),unlist(weapon6[4]),unlist(weapon7[4]),unlist(weapon8[4]))
timeseriesarray <- rbind(timeseriesarray,tobind)

armorhp <- ship[4]
shieldhp <- ship[3]
hullhp <- ship[1]
shieldregen <- ship[2]
shieldmax <- ship[3]
armorhp <- ship[4]
startingarmor <- ship[4]
weapon1shots <- 1
weapon2shots <- 1
weapon3shots <- 1
weapon4shots <- 1
weapon5shots <- 1
weapon6shots <- 1
weapon7shots <- 1
weapon8shots <- 1
armormatrix <- matrix(ship[4]/15,5,ship[6]+4)
startingarmormatrixsum <- sum(armormatrix)
timetokill <- 0
#          }
#        }
#      }
#    }
#  }
#}
# }
}
#}
colnames(timeseriesarray) <-  c("Timetokill", "Weapon1", "Weapon2", "Weapon3", "Weapon4", "Weapon5", "Weapon6", "Weapon7", "Weapon8")

sortbytime <- timeseriesarray[order(as.integer(timeseriesarray$Timetokill)),]

write.table(sortbytime, file = paste("optimizeweaponsbytime",shipname,"allweaponswithacc.txt", sep=""), row.names=FALSE, sep="\t")
}

[close]

This does the italicize things so as a zip in attachment.

Edit to add: there might be an actual error in the running index at the loop at the end as that should probably be (max i + max j + max k ... + p). So if so then that should be fixed. Whether it is so depends on whether I accounted for this in the permutations matrix. If correct then that would have resulted in using improper accuracy parameters and is reason enough to re-compute. Unfortunately simply no time to go through the code and check for now, need to do later. Apologies about a possible error.

E2: yes there is an error. You should replace i with length(weapon1choices) etc from weapon 2 onwards in the running index. This means individual weapon timeseries were correct but the weapon comparison script was assigning incorrect accuracy to weapons much of the time. Well, good thing to have found it. I can probably set my one laptop on re-crunching numbers while doing real research on the other, so expect an update on data and numbers later.

Corrected, this section of the code reads

code

    if(weapon1[4] != ""){
      weapon1[[6]] <- lookup(f,i,1)
      weapon1[[7]] <- lookup(f,i,2)
    }
    
    if(weapon2[4] != ""){
      weapon2[[6]] <- lookup(f,length(weapon1choices)+j,1)
      weapon2[[7]] <- lookup(f,length(weapon1choices)+j,2)
    }
    
    if(weapon3[4] != ""){
      weapon3[[6]] <- lookup(f,length(weapon1choices)+length(weapon2choices)+k,1)
      weapon3[[7]] <- lookup(f,length(weapon1choices)+length(weapon2choices)+k,2)
    }
    
    if(weapon4[4] != ""){
      weapon4[[6]] <- lookup(f,length(weapon1choices)+length(weapon2choices)+length(weapon3choices)+l,1)
      weapon4[[7]] <- lookup(f,length(weapon1choices)+length(weapon2choices)+length(weapon3choices)+l,2)
    }
    
    if(weapon5[4] != ""){
      weapon5[[6]] <- lookup(f,length(weapon1choices)+length(weapon2choices)+length(weapon3choices)+length(weapon4choices)+m,1)
      weapon5[[7]] <- lookup(f,length(weapon1choices)+length(weapon2choices)+length(weapon3choices)+length(weapon4choices)+m,2)
    }
    
    if(weapon6[4] != ""){
      weapon6[[6]] <- lookup(f,length(weapon1choices)+length(weapon2choices)+length(weapon3choices)+length(weapon4choices)+length(weapon5choices)+n,1)
      weapon6[[7]] <- lookup(f,length(weapon1choices)+length(weapon2choices)+length(weapon3choices)+length(weapon4choices)+length(weapon5choices)+n,2)
    }
    if(weapon7[4] != ""){
      weapon7[[6]] <- lookup(f,length(weapon1choices)+length(weapon2choices)+length(weapon3choices)+length(weapon4choices)+length(weapon5choices)+length(weapon6choices)+o,1)
      weapon7[[7]] <- lookup(f,length(weapon1choices)+length(weapon2choices)+length(weapon3choices)+length(weapon4choices)+length(weapon5choices)+length(weapon6choices)+o,2)
    }
    if(weapon8[4] != ""){
      weapon8[[6]] <- lookup(f,length(weapon1choices)+length(weapon2choices)+length(weapon3choices)+length(weapon4choices)+length(weapon5choices)+length(weapon6choices)+length(weapon7choices)+p,1)
      weapon8[[7]] <- lookup(f,length(weapon1choices)+length(weapon2choices)+length(weapon3choices)+length(weapon4choices)+length(weapon5choices)+length(weapon6choices)+length(weapon7choices)+p,2)
    }

[close]

[attachment deleted by admin]

[Liral]

Wow, you are fast and smart. Thanks for the effort and the contribution. It seems like this rewritten script does not actually write the necessary CSVs from text files, so that needs to be added.

It wasn't speed or brains but spending hours copying, pasting, fiddling, deleting, wondering, repeating, and searching the internet for R syntax, the table manipulation parts of which still largely elude me, while thinking, "Gosh, what have I gotten myself into?  I hope I won't have to spend as much time on the next two files!", and feeling as though I were navigating a maze with a tin bucket stuck over my head. 

You're welcome for the contribution, though!  This software could show modders how their weapons balance.  I don't understand what making it "write the necessary CSVs from text files" would entail, so would please you elaborate?

Unfortunately some real life research work is becoming urgent. So despite being a little obsessed with Starsector just now I absolutely must set this aside for a couple of days at least and spend free time working on other things (I have somehow ended up trying to do a PhD while also doing two other degrees and working full time, the upside is there is never time to think about why you would ever want to do such a thing).

I think you have filled your days with academic undertakings and that this project is an extra-curricular one.

Anyway I think you mean you need the txt file (really tsv but I like Windows notepad) data containing the results from running the first script. So here it is. https://filetransfer.io/data-package/aWFSUaVb#link

Yes indeed!  Thanks so much.  I can now continue development!    Wooo!

I can't comment on every step of the code just now but here is a short version and I can go in more depth at another timepoint, also with some extra comments on already commented sections

Thanks, I appreciate that you took the time to write these comments, which will help me understand the code!

[CapnHector]

You're welcome for the contribution, though!  This software could show modders how their weapons balance.  I don't understand what making it "write the necessary CSVs from text files" would entail, so would please you elaborate?

Your script requests a file called e.g. medium_guns.csv but does not write one, whereas weaponsoptimize3.R outputs files in .txt. Probably just WIP. E: nevermind, I see, you meant to include the ones in your post.

[CapnHector]

Alright so as planned above I let this run while I was working on some graphs and glms etc. on my other laptop and it's done. The results were not much changed from the original but now the running index problem is corrected. This also fixes the incorrect width on the Dominator. I'll update post #1 with the results. Here is the raw data: https://filetransfer.io/data-package/C1m2MwU8#link , and attached are fixed scripts. This also fixed the medium missiles problem which was pretty trivial.

[attachment deleted by admin]

[Liral]

Thanks for the code and raw data!  I have almost entirely restructured the optimize script and had an idea while working on it: if only we could directly calculate the probability of hitting a cell, then the code could skip the lengthy step of randomly generating hundreds or thousands of angles and errors and comparing them to the corresponding angle bounds.

Reviewing this code to turn it into a formula, I have noticed a portion I might not understand but, if I do, suspect to be incorrect.

Code

#where does one shot hit within the weapon's arc of fire, in pixel difference from the target? get a random angle in degrees according to a uniform distribution, 
#then consider that the circumference is 2 pi * range pixels, so the hit coordinates in pixels are
shotangle <- function(acc) return(runif(1,-acc/2,acc/2)/360*2*pi*range)
#now add a random positional error to the coordinates of the hit
hitlocation <- function(acc){
  location <- shotangle(acc)
  location <- location + rnorm(1,0,error)
  return(location)
}

I would write the function to generate a random shot angle as

Code

shotAngle <- function(spread) {
    return(runif(1, -spread/2, spread/2) * pi / 180))
}

Meanwhile, vectorizing the code and switching from linear to binary search should noticeably speed it up.

Code

#engagementrange
RANGE <- 1000

#fudge factor
ERROR_SD <- 0.05

#the fudge factor should be a function of range (more error in 
#position at greater range), but not a function of weapon 
#firing angle, and be expressed in terms of pixels
ERROR <- ERROR_SD * RANGE

SHOTS <- 100000

getShotAngles <- function(spread) { 
    return(runif(SHOTS, -spread / 2, spread / 2) * pi / 180) 
}

#now add a random positional error to the coordinates of the hit
getShotErrors <- function() {
    return(rnorm(SHOTS, 0, ERROR))
}

#so which box was hit?
isCellHit <- function(angle, intervalBounds, ship) {
    if (angle < intervalBounds[1]) return(1)
    if (angle > tail(1, intervalBounds)) return(ship[6] + 2)
    left <- 1
    right <- length(intervalBounds)
    while(TRUE) {
        index <- floor((left + right) / 2)
        if (angle <= bounds[index]) right <- index
        else if (angle > bounds[index + 1]) left <- index
        else return(index + 1)
    }
}

#generates the shot distribution per shot
getDistribution <- function(spread) {
    distribution <- vector(mode = "double", length = ship[6] + 2)
    hitLocations <- isCellHit(getShotAngles(spread) + getShotErrors())
    distribution[hitLocations] <- distribution[hitLocations] + 1
    return(distribution / sum(distribution))
}

[CapnHector]

It's great to have a real programmer looking at this. See, I've only learned the bits necessary to do statistical stuff and never give a thought to code clarity or what might be an optimal search algorithm. I'll probably learn stuff that I can even use in real life.

Anyway, you are absolutely correct that the shot angle function does not generate a shot angle.

The shot angle, in radians, is runif(1,-acc/2,acc/2)/360*2*pi. However, on the next line we are adding a random positional error (rather than angle error) to this. So despite my function being called shotangle it actually does what it says in the comments and returns arc length (with sign) from 0 to where the shot hits, to which is then added a random positional error and then compared vs the ship 's and individual armor cells' width to find hit location on ship.

You could also generate the distribution analytically. It is given by Y=X_1+X_2 where Y is shot hit coordinate in pixels, X_1 is the uniform distribution U(-arclength, arclength) and X_2 is the normal distribution N(0,errorSD^2). Then compute P(upperbound > Y > lowerbound) for each cell. We had an argument in the previous thread about what is the correct distribution and ended up simulating.

I'll add that if you are considering generalizing this method to arbitrary ranges and ship widths then should probably correct for curvature. I have not done that as the effect is expected to be negligible at range 1000, but really we would be interested in the coordinates of the projection of the edge of each armor cell to the circle defined by range, so if the range is short and the ship is wide then this would be something that needs to be considered. To avoid this issue, you might instead of the arc length conceptualize that the cells and positional error are not along an arc defined by weapon range, but are instead perpendicular to the firing ship. Then instead of calculating arc length we would use tan(shotangle)*range so the position hit is understood to be the height of the right angle triangle defined by shot angle and range.

Like this (excuse my drawing skills):

[attachment deleted by admin]

[Vanshilar]

I have almost entirely restructured the optimize script and had an idea while working on it: if only we could directly calculate the probability of hitting a cell, then the code could skip the lengthy step of randomly generating hundreds or thousands of angles and errors and comparing them to the corresponding angle bounds.

To me that's an important benefit of CapnHector's "probability wave" approach: that all the damages could be calculated in one step. Coding-wise, this would be vectorizing all sorts of steps. I haven't really looked at the code, but the way I coded it in my Excel spreadsheet was:

1. Make use of the fact that across the armor band, the top and bottom cells have the same values, and the middle 3 rows of cells have the same values. Therefore, the armor cell matrix only needs to be 2 rows. This does mean you have to be careful to do *2 and *3 at the appropriate points though, but means you only need to carry around a 2 x N matrix rather than a 5 x N matrix representing the armor cells.

2. Calculate the equivalent armor rating for each hittable armor cell, which is a 1 x (N - 4) vector, since the weapon can never hit the farther 2 armor cells on either side. That's a straightforward formula to relate the armor matrix to this vector, which is:

1/30 * [ 2 * (0 i+1 i+2 i+3 0) of row 1 of armor matrix + 3 * (i 2*(i+1) 2*(i+2) 2*(i+3) i+4) of row 2 of armor matrix]

or the following matrix dotted with each applicable cell of the matrix of armor cells:

0 1/15 1/15 1/15 0
1/10 1/5 1/5 1/5 1/10

noting that the i's are the indices of the elements, and the inner "2*" refers to multiplying the value of that element by 2. But I'm not sure if there's an easy function to do this in one step, or if you end up having to use a for loop for each element. Then each element in this 1 x (N - 4) vector is set to 0.05 * base armor rating if the value is below that.

3. Input the desired probability distribution, which is also a 1 x (N - 4) vector, and does not necessarily sum to 1. (Any amount less than 1 represents a miss.) This is the part you're talking about here.

4. Calculate the hit damage reduction for each hittable armor cell. That's just a 1 x (N - 4) vector of hit_str/(hit_str + equiv_armor_rating), with a minimum value of 0.15, but again I'm not sure if there's an easy coding way to apply it cell-by-cell in one step.

5. Calculate the damage applied if the shot hit each armor cell. This could be combined with step #4 if desired. If it's a separate step, then it's just a simple scalar (weapon damage) multiplying the hit damage reduction vector.

6. Multiply the damage applied in step #5 with the probability vector in step #3. This is just multiplying two vectors together, element by element, resulting in a new 1 x (N - 4) vector that represents the net damage applied in this round if the shot hit that cell.

7. Expand the net damage vector in step #6 into the damage to be applied to each armor cell. So you're expanding a 1 x (N - 4) vector into a 2 x N matrix by applying a convolution (using the matrix above as the kernel). Not sure if R has a shortcut for this, but since convolution is common in image processing, it might have some function or another to do this as one step.

8. Subtract out the damage matrix from the armor matrix. Not sure how easy it is to apply checking for what happens when no more armor is left; in my case, I used it as a way to count up the hull damage.

Anyway, a lot of these computation steps can probably be vectorized, which would speed up computation greatly. Not sure if some of them might already be done that way in CapnHector's code, but it's worth considering if it's not.

I've attached a condensed version of the Excel spreadsheet I used for my earlier calculations; hope it helps. (The original spreadsheet goes down to 500 steps, and also has a second page with similar code except it hits 1 cell of the cells, generated at random, for each of the 500 steps, but it doesn't fit on the forum.) The input and results areas are to the right, in column AD onward, whereas the calculations are done to the left. The user-input hit probability vector are cells E5:X5. The subsequent steps automatically copy the values in those cells if they're modified. It allows for a hittable cell width of up to 20 cells, which means the armor matrix is 24 cells wide. However there's no harm if you use fewer cells than that; what I have set there is 15.

As a side note, in terms of future expansion, it may be worth:

1. Having the minimum damage reduction value be a user-defined parameter and not hard-coded, since Polarized Armor can reduce it from 0.15 to 0.10.
2. Having the ability to decrease incoming damage to armor by a user-defined percentage (i.e. for Impact Mitigation). Note that the way it works is it actively decreases the incoming shot's damage, i.e. if there's armor, then the incoming shot's damage goes from 400 to 300. (A double whammy since it decreases not just the damage but the hit strength as well, which makes Impact Mitigation really good.)
2. Having the ability to decrease incoming damage to hull by a user-defined percentage (i.e. for Damage Control).
4. Having the ability to decrease the armor hit strength by a user-defined percentage (i.e. for the upcoming Ablative Armor hullmod).

These are not necessary but may be useful for playing around in the future with different builds and different targets.

[attachment deleted by admin]

[Liral]

It's great to have a real programmer looking at this.

I wish!  I'm just a hobbyist. 

See, I've only learned the bits necessary to do statistical stuff and never give a thought to code clarity or what might be an optimal search algorithm. I'll probably learn stuff that I can even use in real life.

Imagining a program as a graph of nodes representing functions and data structures and edges representing function calls reveals what makes code well-formed: every group of statements doing one thing is wrapped inside a function, and functions call other functions, which call other functions, and so on until the 'driver' code at the root of the graph, which should be structured to clearly expose the program to a reader who does not know the contents of the functions.  A node can also contain  functions and data structures associated into an objects, which can act on itself or other objects.  Beware of absolutist ideas of how to arrange these programs, for great holy wars have been waged over them, leaving only long-unread flaming comments and clumsy, scarcely-used languages, but rather make your program run to suit your purpose, readable to make it easy to modify, and finally efficient to save time if it is slow but must be run often.

Anyway, you are absolutely correct that the shot angle function does not generate a shot angle.

Please, please name functions after what they do.  I 

The shot angle, in radians, is runif(1,-acc/2,acc/2)/360*2*pi. However, on the next line we are adding a random positional error (rather than angle error) to this. So despite my function being called shotangle it actually does what it says in the comments and returns arc length (with sign) from 0 to where the shot hits, to which is then added a random positional error and then compared vs the ship 's and individual armor cells' width to find hit location on ship.

Now I understand!  Thanks for explaining.

You could also generate the distribution analytically. It is given by Y=X_1+X_2 where Y is shot hit coordinate in pixels, X_1 is the uniform distribution U(-arclength, arclength) and X_2 is the normal distribution N(0,errorSD^2). Then compute P(upperbound > Y > lowerbound) for each cell. We had an argument in the previous thread about what is the correct distribution and ended up simulating.

Running the numbers both ways and comparing might reveal whether the methods are close enough for the user's purposes.  I speculate that the angle from the front left corner of each cell to the weapon to the front right corner of that cell, divided by the angle from the left edge of the ship to the weapon to the right edge of the ship, would approach the per-cell probabilities of large random samples.

I'll add that if you are considering generalizing this method to arbitrary ranges and ship widths then should probably correct for curvature. I have not done that as the effect is expected to be negligible at range 1000, but really we would be interested in the coordinates of the projection of the edge of each armor cell to the circle defined by range, so if the range is short and the ship is wide then this would be something that needs to be considered. To avoid this issue, you might instead of the arc length conceptualize that the cells and positional error are not along an arc defined by weapon range, but are instead perpendicular to the firing ship. Then instead of calculating arc length we would use tan(shotangle)*range so the position hit is understood to be the height of the right angle triangle defined by shot angle and range.

I profiled r*tan(x) and r*x and found the former to take ten times longer, and |r*tan(x) - r*x| is single-digit for angles less than about 10 degrees, though I have found decreasing the range to increase the interval of angles across which the difference is as small.

Like this (excuse my drawing skills):

Yeah, this makes good sense.

[CapnHector]

Anyway, you are absolutely correct that the shot angle function does not generate a shot angle.

Please, please name functions after what they do.  I 

In my defense that is what the function originally did, but then I changed my mind 

You could also generate the distribution analytically. It is given by Y=X_1+X_2 where Y is shot hit coordinate in pixels, X_1 is the uniform distribution U(-arclength, arclength) and X_2 is the normal distribution N(0,errorSD^2). Then compute P(upperbound > Y > lowerbound) for each cell. We had an argument in the previous thread about what is the correct distribution and ended up simulating.

Running the numbers both ways and comparing might reveal whether the methods are close enough for the user's purposes.  I speculate that the angle from the front left corner of each cell to the weapon to the front right corner of that cell, divided by the angle from the left edge of the ship to the weapon to the right edge of the ship, would approach the per-cell probabilities of large random samples.

What do you mean? Isn't that the exact thing the program does? I was too lazy/busy to calculate what the distribution is but fortunately somebody ahs already done it:
https://www.cambridge.org/core/books/abs/measurement-uncertainty-and-probability/the-sum-of-normal-and-uniform-variates/DECD19AD9E3D68916CC4B04532418F58

We can implement this in R as

G <- function(y) return(y*pnorm(y) + dnorm(y))
fEz <- function(z, a, b) return((1/2/b)*(pnorm(z/a+b/a)-pnorm(z/a-b/a)))
PrEltZ <- function(z, a, b) return(a/2/b*(G(z/a+b/a)-G(z/a-b/a)))

so now a cell's hit probability should be given by PrEltZ(cellupperbound, a, b)-PrEltZ(cellowerbound, a, b) where b is the maximum arc length/right angled triangle height for each shot and a is the standard deviation of the normal error.

See also
matrix <- matrix(data=0,nrow = 1000,ncol=1)
for (i in 1:1000) {
matrix[i,1] <- fEz((i-500)/100,1,4)
}
plot(matrix)


Let's see if this produces the same distribution that the program does. For the simulation we have

range<-1000
error <- 1000*0.05
shotangle <- function(acc) return(runif(1,-acc/2,acc/2)/360*2*pi*range)
#now add a random positional error to the coordinates of the hit
hitlocation <- function(acc){
  location <- shotangle(acc)
  location <- location + rnorm(1,0,error)
  return(location)
}

matrix <- matrix(data=0,nrow = 100000,ncol=1)
for (i in 1:100000) {
  matrix[i,1] <- hitlocation(10)
}
hist(matrix[matrix<200 & matrix > -200])

Result


For the mathematical expression we have

matrix <- matrix(data=0,nrow = 400,ncol=1)
for (i in 1:400) {
  matrix[i,1] <- fEz(i-200,error,10/2/360*2*pi*range)
}
plot(matrix)

Result

[Liral]

What do you mean? Isn't that the exact thing the program does? I was too lazy/busy to calculate what the distribution is but fortunately somebody ahs already done it:
https://www.cambridge.org/core/books/abs/measurement-uncertainty-and-probability/the-sum-of-normal-and-uniform-variates/DECD19AD9E3D68916CC4B04532418F58

I meant that I think a simplified analytical model—divide the angle for each cell by the total angle of the ship—would well-enough approximate the numerical one for the users' purposes. I didn't expect covering the subject to entail a big, thick book: statistics is an enormous field! 

We can implement this in R as

G <- function(y) return(y*pnorm(y) + dnorm(y))
fEz <- function(z, a, b) return(pnorm(z/a+b/a)-pnorm(z/a-b/a))
PrEltZ <- function(z, a, b) return(a/2/b*(G(z/a+b/a)-G(z/a-b/a)))

so now a cell's hit probability should be given by PrEltZ(cellupperbound, a, b)-PrEltZ(cellowerbound, a, b) where b is the maximum arc length/right angled triangle height for each shot and a is the standard deviation of the normal error.

Wooo!  Thanks for finding and implementing the relevant functions—how convenient!    Now I need your help ensuring the code is readable and maintainable: would you please tell me whether you want to keep fEz and then accordingly spell-out the names of the functions and that of the variables y and z and fill the bracketed docstring portions?

#[in a sentence, what does this function do?]
#
#[if that sentence isn't enough, start the paragraph for more here and leave a line]
#
#y - [what is y?]
G <- function(y) return(y * pnorm(y) + dnorm(y))

#[in a sentence, what does this function do?]
#
#[if that sentence isn't enough, start the paragraph for more here and leave a line]
#
#z - [what is z?]
#maximumDeviation - distance a shot may deviate from [what?] and still hit [which?] cell
#standardDeviation - standard deviation of [what?]
fEz <- function(z, maximumDeviation, standardDeviation) { 
    return(pnorm((z + maximumDeviation) / standardDeviation) 
            - pnorm((z - maximumDeviation) / standardDeviation))
}

#[in a sentence, what does this function do?]
#
#[if that sentence isn't enough, start the paragraph for more here and leave a line]
#
#z - [what is z?]
#maximumDeviation - distance a shot may deviate from [what?] and still hit [which?] cell
#standardDeviation - standard deviation of [what?]
PrEltZ <- function(z, maximumDeviation, standardDeviation) {
    return(standardDeviation
            / (2 * maximumDeviation) 
            * (G((z + maximumDeviation) / standardDeviation) 
                - G((z - maximumDeviation) / standardDeviation))
}

[CapnHector]

Let's see. I mean PrEltZ is probably not the user friendliest function name so let's try to make that more reasonable

#Calculate where shots will hit on average. The coordinates of an individual hit are E = E1+E2, where E1 is a uniformly distributed variable representing shot angle and E2 is a randomly distributed variable representing error in the ships' position compared to a situation where the enemy ship is in the middle of the weapon's firing arc and not moving, and the weapon targets the middle.
#
#This helper function is used to define the cumulative distribution function used to calculate where shots hit. y is any real number.
G <- function(y) return(y * pnorm(y) + dnorm(y))

#If we want to show the probability density function of E, it is given by the expression f_E(z)=1/(2b)(Phi(z/a+b/a)-Phi(z/a-b/a)), where the uniform distribution is U(-b,b) and the normal distribution is N(0,a^2) (see Willink 2013, pp. 261)
#This is not used in the code but might be useful for bug testing
#maximumDeviation - the maximum signed arc length of the shot given a maximum firing angle and a range, representing weapon accuracy
#standardDeviation - standard deviation of the normal distribution representing positional error
#z is the shot hit location in pixels and fEz is the probability density function of shot hit location given the average assumptions
#fEz <- function(z, maximumDeviation, standardDeviation) { 
#    return(1/(2*maximumDeviation)*pnorm((z + maximumDeviation) / standardDeviation) 
#            - pnorm((z - maximumDeviation) / standardDeviation))
#}

#This function returns the cumulative distribution function of the distribution E at point z. That is, given we have a shot that hits a random location E, what is the probability that the shot's hit coordinates are less than z? 
#maximumDeviation - the maximum signed arc length of the shot given a maximum firing angle and a range, representing weapon accuracy. Note that 0 represents hitting the weapon's target.
#standardDeviation - standard deviation of the normal distribution representing positional error
probabilityHitLocationLessThanZ <- function(z, maximumDeviation, standardDeviation) {
    return(standardDeviation
            / (2 * maximumDeviation) 
            * (G((z + maximumDeviation) / standardDeviation) 
                - G((z - maximumDeviation) / standardDeviation))
}

This part of the code will probably simply require a little knowledge of statistics and unfortunately not much can be done about that. Also note that I forgot the /2b from fEz when posted initially, added here. E: fixed a typo and added a mention of 0 representing weapon hitting target.

[Liral]

Let's see. I mean PrEltZ is probably not the user friendliest function name so let's try to make that more reasonable

 

This part of the code will probably simply require a little knowledge of statistics and unfortunately not much can be done about that. Also note that I forgot the /2b from fEz when posted initially, added here. E: fixed a typo and added a mention of 0 representing weapon hitting target.

I'll do my best.  Also, long lines of code and long comments are more readable with broken into shorter lines.  Also, what would a descriptive name for G and y be?

#Cumulative distribution of hit location
# 
#The location of a hit is the distance in pixels from the
#ship centerline and equals the sum of a uniformly
#random shot angle and a normally random perpendicular
#distance from the ship centerline
#
#angle - angle between the line from the weapon to the ship
#        and the line along which the weapon aims
#
hitCumulativeDistribution <- function(aimAngle) {
    return(aimAngle * pnorm(aimAngle) + dnorm(aimAngle))
}

#Return hit probability density at an aim angle. 
#
#The probability density function of hit coordinate is 
#f(z) = 1 / (2b) * (Phi((z + b) / a) - Phi((z - b) / a)
#where the uniform distribution is U(-b, b) 
#and the normal distribution is N(0, a^2) 
#(see Willink 2013, pp. 261)
#This function is not used but might help bug testing
#
#maximumDeviation - maximum signed arc length of the shot
#                   given a maximum firing angle and a range, 
#                   representing weapon accuracy
#standardDeviation - standard deviation of the normal distribution 
#                    representing positional error
#hitLocation - shot hit location in pixels
#
hitProbabilityDensity <- function(hitLocation, 
                                  maximumDeviation, 
                                  standardDeviation) 
{ 
    return(1/(2 * maximumDeviation) 
            * pnorm((hitLocation + maximumDeviation) / standardDeviation) 
            - pnorm((hitLocation - maximumDeviation) / standardDeviation))
}

#Returns the cumulative distribution function
#
#The cumulative distribution of the hit probability density is, 
#the probability that the hit coordinate is within some width
#of the ship centerline?
#
#width - maximum acceptable hit coordinate
#maximumDeviation - the maximum signed arc length of the
#                   shot given a maximum firing angle and a 
#                   range, representing weapon accuracy.
#                   0 represents hitting the target.
#standardDeviation - standard deviation of the normal 
#                   distribution representing positional error
#
probabilityHitWithinWidth <- function(width, 
                                      maximumDeviation, 
                                      standardDeviation) 
{
    return(standardDeviation
            / (2 * maximumDeviation) 
            * (hitCumulativeDistribution((width + maximumDeviation) / standardDeviation) 
                - hitCumulativeDistribution((width - maximumDeviation) / standardDeviation))
}

I am still a little confused on two points.  First, although I know that hit location and aim angle are almost interchangeable, I am losing track of whether these functions are talking about one or the other exactly.  And why do the cumulative distribution function comments mention a uniform distribution, which is not called in the code?

[CapnHector]

These functions are talking about hit location as it would not make sense to add a random parameter describing location to an angle. The functions do not call a uniform distribution, since the functions describe the distribution generated by the sum E=E1+E2, where E1~U(-b,b) and E2~N(0,sd^2), and according to Bhattacharjee et al. who did the math and published it once this follows the distribution presented ie. it can be expressed in terms of the standard normal distribution as presented. I suppose you could call G distributionPlusDensity but I don't know if that's really more informative and it is hard to see what G exactly represents when haven't done the math. But calculating the convolution of the two distributions to re-do it is non-trivial and I don't have access to the original article right now.

Anyway, to make sure this works, here is another test. This is what I get from running the original script with 100000 samples:

Spoiler

#SHIP
#dominator, hullhp, shieldregen, shieldmax, startingarmor, widthinpixels, armorcells
ship <- c(14000, 500, 10000, 1500, 220, 12)

#engagementrange
range <- 1000

#weaponaccuracy - this will be made a function of time and weapon later. the accuracy of a hellbore is 10
acc <- 10

#fudge factor
errorsd <- 0.05
#the fudge factor should be a function of range (more error in position at greater range), but not a function of weapon firing angle, and be expressed in terms of pixels
error <- errorsd*range

#where does one shot hit within the weapon's arc of fire, in pixel difference from the target? get a random angle in degrees according to a uniform distribution,
#then consider that the circumference is 2 pi * range pixels, so the hit coordinates in pixels are
shotangle <- function(acc) return(runif(1,-acc/2,acc/2)/360*2*pi*range)

#how much is the visual arc of the ship in rad?
shipangle <- ship[5]/(2* pi *range)

#how much is the visual arc of a single cell of armor in rad?
cellangle <- shipangle/ship[6]

#now assume the weapon is targeting the center of the ship's visual arc and that the ship is in the center of the weapon's firing arc
#which cell will the shot hit, or will it miss?
#call the cells (MISS, cell1, cell2, ... ,celli, MISS) and get a vector giving the (maximum for negative / minimum for positive) angles for hitting each
anglerangevector <- vector(mode="double", length = ship[6]+1)
anglerangevector[1] <- -shipangle/2
for (i in 1:(length(anglerangevector)-1)) anglerangevector[i+1] <- anglerangevector+cellangle

#now convert it to pixels
anglerangevector <- anglerangevector*2*pi*range

#this vector will store the hits
shipcellvector <- vector(mode="double", length = ship[6]+2)

#now add a random positional error to the coordinates of the hit
hitlocation <- function(acc){
  location <- shotangle(acc)
  location <- location + rnorm(1,0,error)
  return(location)
}

#so which box was hit?
cellhit <- function(angle){
  if(angle < anglerangevector[1]) return(1)
  if(angle > anglerangevector[ship[6]+1]) return(ship[6]+2)
  for (i in 1:length(anglerangevector)) {
    if ((angle > anglerangevector) & (angle <= anglerangevector[i+1])) return(i+1)
  }
}


# this function generates the shot distribution per 1 shot with 100000 samples
createdistribution <- function(acc){
  distributionvector <- vector(mode="double", length = ship[6]+2)
  for (i in 1:100000){
    wherehit <- cellhit(hitlocation(acc))
    distributionvector[wherehit] <- distributionvector[wherehit] +1
  }
  return(distributionvector/sum(distributionvector))
}

# this is the default distribution of damage to armor cells
b <- matrix(0,nrow=5,ncol=5)
b[1:5,2:4] <- 1/30
b[2:4,1:5] <- 1/30
b[2:4,2:4] <- 1/15
b[1,1] <- 0
b[1,5] <- 0
b[5,1] <- 0
b[5,5] <- 0

#this function generates a sum of matrices multiplied by the distribution

createhitmatrix <- function(acc){
  hitmatrix <- matrix(0,5,ship[6]+4)
  distributionvector <- createdistribution(acc)
  for (i in 1:ship[6]){
    hitmatrix[,i:(i+4)] <- hitmatrix[,i:(i+4)]+b*(distributionvector[i+1])
  }
  return(hitmatrix)
}


plot(createdistribution(10))

[close]

and here is what I get from running this:

probabilityvector <- vector(mode="double",14)
for (i in 2:length(probabilityvector)){
    probabilityvector[\i] <- PrEltZ(anglerangevector[\i],error,10/2/360*2*pi*range)-PrEltZ(anglerangevector[i-1],error,10/2/360*2*pi*range)
}
probabilityvector[1] <- PrEltZ(anglerangevector[1],error,10/2/360*2*pi*range)
probabilityvector[14] <- 1-PrEltZ(anglerangevector[13],error,10/2/360*2*pi*range)
plot(probabilityvector)

(note, escaped i due to bbcode)


But the mathematical method is very much faster.

[Liral]

These functions are talking about hit location as it would not make sense to add a random parameter describing location to an angle.

Ok, I have changed angle references to position references.

The functions do not call a uniform distribution, since the functions describe the distribution generated by the sum E=E1+E2, where E1~U(-b,b) and E2~N(0,sd^2), and according to Bhattacharjee et al. who did the math and published it once this follows the distribution presented ie. it can be expressed in terms of the standard normal distribution as presented.

I have changed the documentation and structure of the third function accordingly but wonder if it is clear and accurate.  Please tell me what changes I should make further.

I suppose you could call G distributionPlusDensity but I don't know if that's really more informative and it is hard to see what G exactly represents when haven't done the math. But calculating the convolution of the two distributions to re-do it is non-trivial and I don't have access to the original article right now.

I have inlined the first function into the third.  It is ugly but perhaps less-so than the helper function.

Anyway, to make sure this works, here is another test. This is what I get from running the original script with 100000 samples:

Spoiler

#SHIP
#dominator, hullhp, shieldregen, shieldmax, startingarmor, widthinpixels, armorcells
ship <- c(14000, 500, 10000, 1500, 220, 12)

#engagementrange
range <- 1000

#weaponaccuracy - this will be made a function of time and weapon later. the accuracy of a hellbore is 10
acc <- 10

#fudge factor
errorsd <- 0.05
#the fudge factor should be a function of range (more error in position at greater range), but not a function of weapon firing angle, and be expressed in terms of pixels
error <- errorsd*range

#where does one shot hit within the weapon's arc of fire, in pixel difference from the target? get a random angle in degrees according to a uniform distribution,
#then consider that the circumference is 2 pi * range pixels, so the hit coordinates in pixels are
shotangle <- function(acc) return(runif(1,-acc/2,acc/2)/360*2*pi*range)

#how much is the visual arc of the ship in rad?
shipangle <- ship[5]/(2* pi *range)

#how much is the visual arc of a single cell of armor in rad?
cellangle <- shipangle/ship[6]

#now assume the weapon is targeting the center of the ship's visual arc and that the ship is in the center of the weapon's firing arc
#which cell will the shot hit, or will it miss?
#call the cells (MISS, cell1, cell2, ... ,celli, MISS) and get a vector giving the (maximum for negative / minimum for positive) angles for hitting each
anglerangevector <- vector(mode="double", length = ship[6]+1)
anglerangevector[1] <- -shipangle/2
for (i in 1:(length(anglerangevector)-1)) anglerangevector[i+1] <- anglerangevector+cellangle

#now convert it to pixels
anglerangevector <- anglerangevector*2*pi*range

#this vector will store the hits
shipcellvector <- vector(mode="double", length = ship[6]+2)

#now add a random positional error to the coordinates of the hit
hitlocation <- function(acc){
  location <- shotangle(acc)
  location <- location + rnorm(1,0,error)
  return(location)
}

#so which box was hit?
cellhit <- function(angle){
  if(angle < anglerangevector[1]) return(1)
  if(angle > anglerangevector[ship[6]+1]) return(ship[6]+2)
  for (i in 1:length(anglerangevector)) {
    if ((angle > anglerangevector) & (angle <= anglerangevector[i+1])) return(i+1)
  }
}


# this function generates the shot distribution per 1 shot with 100000 samples
createdistribution <- function(acc){
  distributionvector <- vector(mode="double", length = ship[6]+2)
  for (i in 1:100000){
    wherehit <- cellhit(hitlocation(acc))
    distributionvector[wherehit] <- distributionvector[wherehit] +1
  }
  return(distributionvector/sum(distributionvector))
}

# this is the default distribution of damage to armor cells
b <- matrix(0,nrow=5,ncol=5)
b[1:5,2:4] <- 1/30
b[2:4,1:5] <- 1/30
b[2:4,2:4] <- 1/15
b[1,1] <- 0
b[1,5] <- 0
b[5,1] <- 0
b[5,5] <- 0

#this function generates a sum of matrices multiplied by the distribution

createhitmatrix <- function(acc){
  hitmatrix <- matrix(0,5,ship[6]+4)
  distributionvector <- createdistribution(acc)
  for (i in 1:ship[6]){
    hitmatrix[,i:(i+4)] <- hitmatrix[,i:(i+4)]+b*(distributionvector[i+1])
  }
  return(hitmatrix)
}


plot(createdistribution(10))

[close]

and here is what I get from running this:

probabilityvector <- vector(mode="double",14)
for (i in 2:length(probabilityvector)){
    probabilityvector[\i] <- PrEltZ(anglerangevector[\i],error,10/2/360*2*pi*range)-PrEltZ(anglerangevector[i-1],error,10/2/360*2*pi*range)
}
probabilityvector[1] <- PrEltZ(anglerangevector[1],error,10/2/360*2*pi*range)
probabilityvector[14] <- 1-PrEltZ(anglerangevector[13],error,10/2/360*2*pi*range)
plot(probabilityvector)

(note, escaped i due to bbcode)


But the mathematical method is very much faster.

Fantastic!  Finding this an analytical method spares waiting for the code to finish or optimizing it.  Here's the code redone.  Thanks for helping me make it readable! 

Spoiler

#Return hit probability density for a distance
#from the ship centerline.
#
#The hit probability density a at distance c
#the ship centerline is the sum of the uniform 
#distribution U(-b, b) and normal distribution
#N(0, a^2), and this sum is equivalent to
#
#(Phi((c + b) / a) - Phi((c - b) / a)) / 2b
#
#where Phi is the cumulative density function
#(see Willink 2013, pp. 261).
#
#Note: This function is not used but might help 
#bug-testing.
#
#position - in pixels from centerline of ship
#spread - maximum signed arc length of the shot,
#         given maximum firing angle and range, 
#         representing weapon accuracy
#standardDeviation - standard deviation of the  
#                    normal distribution,  
#                    representing positional 
#                    error
#
hitProbabilityDensity <- function(position, 
                                  spread, 
                                  standardDeviation) 
{ 
    return(1/(2 * spread) 
            * pnorm((position + spread) / standardDeviation) 
            - pnorm((position - spread) / standardDeviation))
}


#Return cumulative distribution of hit 
#probability density.
#
#The cumulative distribution of hit probability
#density is the probability the hit coordinate
#is within some number of pixels of the ship 
#centerline.
#
#width - maximum acceptable hit coordinate
#spread - maximum signed arc length of the
#         shot, given maximum firing angle and 
#         range, representing weapon inaccuracy
#standardDeviation - standard deviation of the 
#                    normal distribution, 
#                    representing positional 
#                    error
#
hitProbabilityWithinWidth <- function(width, 
                                      spread, 
                                      standardDeviation) 
{
    left <- (width + spread) / standardDeviation
    right <- (width - spread) / standardDeviation
    left <- (left * pnorm(left) + dnorm(left)
    right <- (right * pnorm(right) + dnorm(right))
    return(standardDeviation / (2 * spread) * (right - left)
}

[close]