#--Redshift of M77/NGC 1068:
z <- 0.003793

#--NASA Extragalactic Database (NED) URL:
url <- "http://nedwww.ipac.caltech.edu/spc1/1995ApJS...98..129K"

#--Names of ASCII-format optical spectra of NGC1068:
files <- c("NGC_1068:S:B:ksv1995_NED.txt", "NGC_1068:S:R:ksv1995_NED.txt")

#--Define known spectral lines:
lines.lambda <- c(4358, 4363, 4686, 4861, 4959, 5007, 6301, 6548, 6562.81, 6723)
lines.name <- c("Hg iII", "OIII", "He II",
                as.expression(expression(paste("H", beta))), "OIII", "OIII", "OI",
                "NII", as.expression(expression(paste("H", alpha))), "SII")

#--Specify X & Y offsets for line labels w.r.t. line wavelength:
adj.x.lines <- c(0.5, 0.5, 0.5, 0.5, 1, -0.2, 0.5, 1.3, -0.6, 0.5)
adj.y.lines <- c(-0.5, -2, -1, -0.5, -4, -1, -0.5, -1, -1, -1)

#--Adjust lines to galaxy redshift:
lines.lambda <- (1+z)*lines.lambda

#--Speed of light in a vacuum:
c_mps <- 2.9979245620E8        # [m/s]

#--Axis labels:
xlab <- as.expression( expression(paste("Wavelength (", ring(A), ")", sep="") ) )
ylab <- as.expression(expression(paste("Intensity (", 10^{-28}, " ", W/m^2/Hz, ")")))

#-----Loop over spectral data files, read in data and concatenate together:
spec <- NULL		# object (will become data frame) exists but with zero contents
i <- 0
for ( f in files ) {
  i <- i + 1
  #--Add URL path to file name:
  f <- paste(url, f, sep="/")
  #--Read in spectral data (in fixed width format):
  a <- read.fwf(f, widths=c(73, 14, 6, 14))
  a <- data.frame(lambda=c_mps/a$V2/1e-10, intensity=a$V4/1e-28) # Extract relevant columns
  spec <- rbind(spec, a)	# Concatenate data frame: bind rows onto existing data frame
}

#--Change output device to pdf file:
pdf(file="M77_opt_spectrum.pdf", height=4, width=6.5, pointsize=16)

#---Remove extra top margin:
par(mar=c(2.8, 2.8, 0.5, 0.5))  # Trim margin around plot [bottom, left, top, right]
mgp <- c(0.8, 0, 0)  # default c(3,1,0) [title,labels,line]
tcl <- 0.5
par(xaxs="r", yaxs="r")  # Extend axis limits by 4% ("i" does no extension)
par(lab=c(10, 10, 7))
lwd <- 2.5  # Define global line width

#--Set up plot window:
plot(spec, type="l", xlab=xlab, ylab=ylab, main=NULL, tcl=tcl, mgp=mgp,
     cex.axis=0.7, cex.lab=0.7)

#---Axis tickmarks & labels:
majat <- seq(4000, 8000, 1000)		# Specify major X axis tick marks (i.e. numbered ticks)
minat <- seq(3500, 8500, 100)		# Specify minor X axis tick marks
majat.s <- seq(0, 50, 10)		# Major side (Y) axis ticks
minat.s <- seq(0, 50, 1)		# Minor Y axis ticks
par(tcl=0.22)				# Tick length
axis(1, at=majat, labels=F)		# Bottom axis major axis tick marks
axis(1, at=minat, labels=F)		#   "     "   minor axis tick marks
axis(2, at=majat.s, labels=F)		# Left axis major axis tick marks
axis(2, at=minat.s, labels=F)		#   "   "   minor axis tick marks
axis(3, at=minat, labels=F)		# Top axis major axis tick marks
axis(3, at=majat, labels=F, tcl=tcl)	#  "   "   minor axis tick marks
axis(4, at=minat.s, labels=F)		# Right axis major axis tick marks
axis(4, at=majat.s, labels=F, tcl=tcl)	#   "    "   minor axis tick marks

#--Identify lines with "rug" (short lines extending from axis):
rug(lines.lambda, side=3, col="blue", lty=2, ticksize=0.1, lwd=1.5)

#--Loop over lines & annotate plot with name label:
for ( i in 1:length(lines.lambda) ) {
  x <- lines.lambda[i]
  y <- spec$intensity[which.min(abs(spec$lambda-x))]
  text(x, y, lines.name[i], cex=0.5, col="red", adj=c(adj.x.lines[i], adj.y.lines[i]))
}

#--Ensure graphics device finishes cleanly:
dev.off()