# aug 26

setwd("/Users/jameshanley/Documents/work/Espionage")

source("ExtractLinesFromLargeFileFunction.R")

result = ExtractLinesFromLargeFileFunction("ATLASfig1.ps",90000,10000)

LENGTH=unlist(result[3]) ; LENGTH ; summary(LENGTH)

n.l = length(LENGTH) ; n.l

X = result[1][[1]] ; Y = result[2][[1]] ; length(X)

x.min=unlist(result[4]) ; y.min=unlist(result[5])
x.max=unlist(result[6]) ; y.max=unlist(result[7])
c(x.min,y.min,x.max,y.max)

x.range=c(400,530) ; y.range = c(609,685)

plot(x.range,y.range,col="white", xlim=x.range, ylim=y.range)
for(i in 1:length(X) )  lines(X[[i]],Y[[i]] )


# focus on tick marks 
plot(x.range,y.range,col="white", xlim=x.range, ylim=y.range)
for(i in 1:47 ) {
	L = LENGTH[i]	
	xa=0; if (X[[i]][1] > X[[i]][L]) xa=1
	ya=0; if (Y[[i]][1] > Y[[i]][L]) ya=1
	if(L==2) text(X[[i]][1],Y[[i]][1], toString(i), cex=0.5,adj=c(xa,ya) ) 
	if(L==2) lines(X[[i]],Y[[i]]) ;
	j = 1 + (i - 5*floor(i/5))
	lines(X[[i]],Y[[i]], lwd=0.5) 
    if(L==2) text(X[[i]][L],Y[[i]][L], toString(i), cex=0.5,adj=c(xa,ya) ) 
}

i.x.720 = 45; i.x.0 = 37
x1 = ( X[[i.x.720]][1] - X[[i.x.0]][1] ) / 720 ; x1

i.y.12 = 30; i.y.0=36
y1 =  (Y[[i.y.12]][1] - Y[[i.y.0]][1])/0.12  ; y1


ci.0 = (Y[[48]] - Y[[i.y.0]][1])/y1 ; ci.0 ; length(ci.0)
t.0  = (X[[48]] - X[[i.x.0]])/x1 ; t.0 ;

plot(t.0,ci.0,type="l",col="grey70")

ci.1 = (Y[[49]] - Y[[i.y.0]][1])/y1 ; ci.1 ; length(ci.1)
t.1 = (X[[49]] - X[[i.x.0]])/x1 ; t.1
lines(t.1,ci.1,type="l")

s.1=1-ci.1
ratio.1=s.1[2:721]/s.1[1:720]
round(1/(1-ratio.1),0)

every=90
T=seq(every+1,721,every)-every/2
CI.1=ci.1[T] ; i.1 = CI.1 - c(0,CI.1[1:(length(CI.1)-1)]) ; i.1 ; sum(i.1)
CI.0=ci.0[T] ; i.0 = CI.0 - c(0,CI.0[1:(length(CI.0)-1)]) ; i.0 ; sum(i.0)
round(i.1/i.0,2)

n.1=approx(seq(1,721,90),c(10229,8817,7797,6324,5137,3967,2830,1747,831), T)$y
n.0=approx(seq(1,721,90),c(5113,4437,3974,3253,2664,2059,1460,878,421), T)$y

c.1 = n.1*i.1 ; c.1 ; sum(c.1); 
c.0 = n.0*i.0; c.0 ; sum(c.0)

rr=i.1/i.0
mme = exp(1.96*sqrt(1/c.1+1/c.0))
plot(T/365,   rr, ylim=c(0,2), xlim=c(0,2) )
segments(0, 1, 720,1 )
for (i in 1:length(i.1)) {
	segments(T[i]/365, rr[i]/mme[i], T[i]/365, rr[i]*mme[i] ) 
	text(T[i]/365, rr[i], toString(round(c.1[i]+c.0[i],0)), adj=c(-0.25,0.5) ) 
	}
rr
round(100*(1-rr),0)

sum(c.1)/sum(c.0)

X.sd = as.vector(sapply(X,function(x) sd(x)))
Y.sd = as.vector(sapply(Y,function(x) sd(x)))

0.5*sum(c.1[2:8]) / sum(c.0[2:8])

0.5*sum(c.1) / sum(c.0)