psname=as.vector(Sys.getenv("ROUTPUT"))
postscript(file=psname, horizontal=FALSE)
par(mfrow=c(2,1))
par(cex=0.7)
par(las=1)
par(lab=c(5,10,0))
f=data.frame(
type=c(
"HSDL4230",
"HSDL4220",
"HPWT-BD00-F4000",
"HPWT-BD00-E4000",
"HPWT-BH00-H4000",
"HPWT-BH00-G4000",
"HPWT-BD00-D4000",
"HPWT-BH00-F4000",
"HPWT-BD00(-0000)",
"HPWT-BH00-E4000",
"HPWT-BH00-D4000",
"HPWT-BH00(-0000)"),

altname=c(
"HSDL\n4230",
"HSDL\n4220",
"HPWT\n-BD00\n-F4000",
"HPWT\n-BD00\n-E4000",
"HPWT\n-BH00\n-H4000",
"HPWT\n-BH00\n-G4000",
"HPWT\n-BD00\n-D4000",
"HPWT\n-BH00\n-F4000",
"HPWT\n-BD00(\n-0000)",
"HPWT\n-BH00\n-E4000",
"HPWT\n-BH00\n-D4000",
"HPWT\n-BH00(\n-0000)"),

count=c(34.32,19.36,15.1,12.6,11.7,10.2,10.1,8.76,7.55,7.3,5.84,4.38))
f$type=as.character(f$type)
f$altname=as.character(f$altname)

f$gain=10*log10(f$count/1000) # Input is in Gq, we want Tq

barplot(f$count, legend.text=f$type, xpd=FALSE, cex.names=0.4,
names.arg=f$type, col=rainbow(length(f$type)), main="Fig. 1: LED on-axis photon flux intensity ", xlab="LED type",
ylab="Intensity [Gq/sym/sr]", ylim=c(0,36), tck=1, )

barplot(f$gain, legend.text=f$type, cex.names=0.4, xpd=FALSE,
names.arg=f$type, col=rainbow(length(f$type)), main="Fig. 2: LED on-axis photon flux intensity [dBTq/sym/sr]",  xlab="LED type",
ylab="Intensity [dBTq/sym/sr]",ylim=c(-24,-14), tck=1, )

# write.table(f, file="bla")

g<-data.frame(
type=c(60,75,80,90,100,130),
gain=c(29.8,29.52,31.91,31.38,32.28,33.74),
g4220=c(25.56,22.51,27.30,27.2,28.81,29.41),
g4230=c(22.82,19.54,24.02,27.24,25.91,27.21)
)
cols=rainbow(length(g$type))
barplot(g$gain, legend.text=g$type, xpd=FALSE,
names.arg=g$gain, main="Fig. 3: Lens gain with HPWT-B type LEDs", sub="For visible LEDs only (red colour), in dB, loupe lenses", xlab="Loupe diameter [mm]",
ylab="Gain [dB]", ylim=c(29,34), tck=1, col=cols)

barplot(g$g4220, legend.text=g$type, xpd=FALSE,
names.arg=g$g4220, main="Fig. 4: Lens gain with HSDL4220 IR LED [dB]", xlab="Loupe diameter [mm]",
ylab="Gain [dB]", tck=1, ylim=c(22,30), col=cols)

barplot(g$g4230, legend.text=g$type, xpd=FALSE,
names.arg=g$g4230, main="Fig. 5: Lens gain with HSDL4230 IR LED [dB]", xlab="Loupe diameter [mm]",
ylab="Gain [dB]", tck=1, ylim=c(18,28), col=cols)

to=outer(f$gain,g$gain,"+")
to[1,]<-g$g4230+f$gain[1]
to[2,]<-g$g4220+f$gain[2]

fun<-function(th, ph, fig){
persp(,,t(to), theta=th, phi=ph, ltheta=135+90+60, col=rgb(1,0.15,0), shade=0.5,
xlab="Lens (60 to 130)", ylab="Diode (strong to weak)",
zlab="Launch int. [dB]", main=paste("Fig. ",fig,": Launch intensity vs. diode and lens type"),
 cex.axis=0.6)
 }
fun(300,45,6)
fun(155,35,7)
fun(30,35,8)

fun=function(to,flag,fig,sub1){
par(xaxt="n")
par(yaxt="n")
image(,,to, xlab="Diode", ylab="Lens", main=paste("Fig. ",fig,": Launch intensity of lens+diode combination [dB]"), sub=sub1, col=rainbow(200, s=0.8, end=0.5))


# horizontal numbers
len=length(f$type)
for (d in 1:len)
text((d-1)/(len-1),-0.04,f$altname[d])

# vertical numbers
len=length(g$type)
for (d in 1:len)
text(-0.03,(d-1)/(len-1),g$type[d])

if (flag == 1)
contour(,,to, add=TRUE)
else{
	xlen=length(f$type)
	ylen=length(g$type)
	for (d in 1:xlen){
		for (e in 1:ylen){
		text((d-1)/(xlen-1)+0.022,0.05+(e-1)/(ylen-1),round(to[d,e],2))
		}
	}
}
}
fun(to,1,9,"Single transmitter")
fun(to,0,10,"Single transmitter")
fun(to+10*log10(2),0,11,"Dual transmitter")
fun(to+10*log10(3),0,12,"Triple transmitter")

rec<-data.frame(type=c("HPWT-BD00-F4000, 130mm", "HSDL4230, 90mm","HSDL4220, 80mm", "HSDL4220, 60mm"),
gain=c(to[3,6], to[1,4], to[2,3], to[2,1]))
rec$type=as.character(rec$type)

par(xaxt="s")
par(yaxt="s")

barplot(to[1,], legend.text=g$type, xpd=FALSE,
names.arg=g$type, main="Fig. 13: Launch intensity with HSDL4230 IR LED", xlab="Lens type",
ylab="Launch intensity [dBTq/sym/sr]", ylim=c(4,13), tck=1,
col=heat.colors(length(g$type)))

barplot(to[2,], legend.text=g$type, xpd=FALSE,
names.arg=g$type, main="Fig. 14: Launch intensity with HSDL4220 IR LED", xlab="Lens type",
ylab="Launch intensity [dBTq/sym/sr]", ylim=c(4,13), tck=1,
col=heat.colors(length(g$type)))

barplot(to[3,], legend.text=g$type, xpd=FALSE,
names.arg=g$type, main="Fig. 15: Launch intensity with HPWT-BD00-F4000 LED", xlab="Lens type",
ylab="Launch intensity [dBTq/sym/sr]", ylim=c(11,16), tck=1,
col=heat.colors(length(g$type)))

barplot(to[4,], legend.text=g$type, xpd=FALSE,
names.arg=g$type, main="Fig. 16: Launch intensity with HPWT-BD00-E4000 LED", xlab="Lens type",
ylab="Launch intensity [dBTq/sym/sr]", ylim=c(10,15), tck=1,
col=heat.colors(length(g$type)))

barplot(rec$gain, legend.text=rec$type, xpd=FALSE,
names.arg=rec$type, main="Fig. 17: Launch intensity of recommended IR diode + loupe combos", sub="HPWT-BD00-F4000 and 130mm loupe is for comparison", xlab="LED type",
ylab="Launch intensity [dBTq/sym/sr]", ylim=c(8,16), tck=1,
col=topo.colors(length(rec$type)))