# These are the R commands for the analysis of the Richards, et al.
# data set, using the method described in Schafer...

# The data files should be downloaded directly from the ApJ web site,
#
#    www.journals.uchicago.edu/AJ/journal/issues/v131n6/205010/205010.html
#
# and labelled "RichardsTable5.dat", etc. Tables 4, 5, and 6 are required.
# Do not alter the files (don't remove the headers)

# The easiest way to run this is in "batch" mode via the command
#
# R CMD BATCH DR3Commands
#
# at the shell. The result will be an object named results_[lambdax]_[lambday]
#
# R can then be started, and this result file can be loaded via, e.g.,
# > load("results_0.05_0.17")
# The object will be named "results"


source("BivTruncRcomm")

# -----------------------------------------------------------------------------------
# Set the main parameters

# These define the boundaries of the truncation region

threecut = 3.0
imagmax = 20.2
imagmin = 15.0
minredsh = 0.1
maxredsh = 5.30
minabsmag = -30.7
maxabsmag = -23.075

# The list of (lambdax,lambday) combinations to try
# Note that every possible combination will be fit

lambdax = c(0.05)
lambday = c(0.17)


# If the analysis has already completed, you can skip to producing the
# plots, by setting plotsonly = TRUE

plotsonly = FALSE


# Other parameters

reductfactor = 5
res2 = 500
deg = 1


tran = function(x)
{
   log10(pmax(x,0))
}


res = res2/reductfactor


# -----------------------------------------------------------------------------------
# Read in data files

quasdat = read.table("RichardsTable5.dat",skip=23)
Kcorr = read.table("RichardsTable4.dat",skip=22)
RichLumfunc = read.table("RichardsTable6.dat",skip=31)

redsh = quasdat[,3]
mag = quasdat[,4]
absmag = quasdat[,5]



# -----------------------------------------------------------------------------------
# Create the res2 by res2 grid to be used in the analysis

redshgridres = minredsh + 0:(res2-1)*(maxredsh-minredsh)/res2 + 
   (maxredsh-minredsh)/res2/2
absmaggridres = minabsmag + 0:(res2-1)*(maxabsmag-minabsmag)/res2 + 
   (maxabsmag-minabsmag)/res2/2


# -----------------------------------------------------------------------------------
# Find K-Correction for each observation

redsh = redsh + 0.0000001
ObsKcorr = (ceiling(redsh*100)-redsh*100)*Kcorr[floor(redsh*100)+1,2]+
           (redsh*100-floor(redsh*100))*Kcorr[ceiling(redsh*100)+1,2]
redsh = redsh - 0.0000001


dVc = function(z,omegam=0.3,omegalam=0.7,hubble=70)
{
   sol = 2.99792458*(10^5)
   DM = approxint(z,omegam,omegalam,hubble)
   sol/hubble*(DM^2)/sqrt(omegam*((1+z)^3)+omegalam)
}


# -----------------------------------------------------------------------------------
# Function to approximate the integral needed in the conversion from apparent to
# absolute magnitudes
# This gives D_M

approxint = function(z,omegam=0.3,omegalam=0.7,hubble=70)
{
   z = z + 0.00000001
   sol = 2.99792458*(10^5)
   zseq = seq(0,max(z),length=1000)
   hold = 1/(omegam*((1+zseq)^3)+omegalam)^0.5
   for(i in 2:length(hold))
   {
      hold[i] = hold[i] + hold[i-1]
   }
   hold2 = hold
   hold2 = rank(c(z,zseq))[1:length(z)] -rank(z)
   hold[hold2]*max(z)/length(hold)/hubble*sol
}

# This function converts apparent magnitudes to absolute magnitudes

mtoM = function(m,redsh,Kcorr)
{
    m - 5*log10(approxint(redsh)*(1+redsh)/10*1000000) - Kcorr
}

Mtom = function(M,redsh,Kcorr)
{
    M + 5*log10(approxint(redsh)*(1+redsh)/10*1000000) + Kcorr
}


# -----------------------------------------------------------------------------------
# Subset the data set

keep = redsh < maxredsh & redsh > minredsh & absmag < maxabsmag & absmag > minabsmag
keep = keep & (mag <= 19.1 | redsh > threecut) & mag <= imagmax


# -----------------------------------------------------------------------------------
# Create the matrix mask

# Find K-Correction for each redshift grid value

wt1 = ceiling(redshgridres*100)-redshgridres*100
gridKcorr = (wt1)*Kcorr[floor(redshgridres*100)+1,2]+(1-wt1)*Kcorr[ceiling(redshgridres*100)+1,2]

bnd191 = mtoM(19.1,redshgridres,gridKcorr)
bnd150 = mtoM(15.0,redshgridres,gridKcorr)
bnd202 = mtoM(20.2,redshgridres,gridKcorr)
expgrid = expand.grid(redshgridres,absmaggridres)

mask = matrix(T,nrow=res2,ncol=res2)
for(i in 1:res2)
{
   if(redshgridres[i] < threecut)
   {
      mask[i,] = absmaggridres < bnd191[i] & absmaggridres > bnd150[i]
   } else
   {
      mask[i,] = absmaggridres < bnd202[i] & absmaggridres > bnd150[i]
   }
}


# -----------------------------------------------------------------------------------
# Call the function BivTrunc()

X = redsh[keep]
Y = absmag[keep]
xlim = c(minredsh,maxredsh)
ylim = c(minabsmag,maxabsmag)

datwghts = pmin(1/quasdat[keep,7],3)

if(!plotsonly)
{
   results = BivTrunc(X,Y, xlim, ylim, mask,lambdax=lambdax,
        lambday=lambday,verbose=T,deg=deg,datwght=datwghts,
        tau=0,reductfactor=reductfactor)
   save(results,file=paste("results_",lambdax,"_",lambday,sep=""))
}



load(file="results_0.05_0.17")

# Plot the bivariate estimate

levs = ((0:20) + 0.4)/20

contlwd = 1

postscript(file="bivest.eps",width=8,height=8,horiz=F)
contcol = rgb(0,0,1,max=1.0)
contcol = rgb(0,0,0)
contcol2 = rgb(0,0,1,max=1.0)
contcol2 = rgb(0,0,0)
bndcol = "black"
lo1 = round(minredsh*100) + 1
hi1 = round(threecut*100) + 1
hi2 = round(maxredsh*100) + 1
plot(redsh[keep],absmag[keep],pch=".",col=rgb(0.7,0.7,0.7),xlab="Redshift",ylab="Absolute Magnitude",
    ylim=c(maxabsmag,minabsmag),xlim=c(minredsh,maxredsh),cex.lab=1.2,cex.axis=1.2)
lines(Kcorr[lo1:hi1,1],pmin(mtoM(19.1,Kcorr[lo1:hi1,1],Kcorr[lo1:hi1,2]),
    maxabsmag),lwd=3,type="l",lty=2,col=bndcol)
lines(Kcorr[lo1:hi2,1],pmax(pmin(mtoM(15.0,Kcorr[lo1:hi2,1],Kcorr[lo1:hi2,2]),
    maxabsmag),minabsmag),lwd=3,type="l",lty=2,col=bndcol)
lines(Kcorr[hi1:hi2,1],mtoM(20.2,Kcorr[hi1:hi2,1],Kcorr[hi1:hi2,2]),lwd=3,type="l",
    lty=2,col=bndcol)
lines(c(threecut,threecut),c(mtoM(20.2,Kcorr[hi1,1],Kcorr[hi1,2]),
    mtoM(19.1,Kcorr[hi1,1],Kcorr[hi1,2])),lwd=3,type="l",lty=2,col=bndcol)
lines(c(maxredsh,maxredsh),c(mtoM(20.2,Kcorr[hi2,1],Kcorr[hi2,2]),
    min(mtoM(19.1,Kcorr[hi2,1],Kcorr[hi2,2]),minabsmag)),lwd=3,type="l",lty=2,col=bndcol)
contour(results$xseqbivest,results$yseqbivest,results$bivest,add=T,levels=levs,lwd=contlwd,col=contcol)
dev.off()


# Plot the redshift marginal

subseq = seq(3,res,by=5)

ztoplot = redshgridres[subseq]
shft = log10(sum(datwghts))
shft2 = shft - log10(dVc(results$xseqcov)*1622/((180/pi)^2))
shft = shft - log10(dVc(results$xseqbivest)*1622/((180/pi)^2))

postscript(file="redshift_marginal.eps",width=8,height=8,horiz=F)
plotmarginal(results,1,varname="Redshift",trans=log10,cex=1.2,shft=shft,shft2=shft2,linecol="black",errcol="black",
      ylabel=expression(log(rho(z)) (M < -23.075) (Mpc^{-3})))
dev.off()


# Make the plot with all luminosity functions together

condplot = unique(RichLumfunc[,1])
hght = c(rep(1,9),rep(1.5,3))
widt = c(1.1,1,1,1.1,1,1,1.1,1,1,1.1,1)
widt = c(1.2,1,1)
hght = c(1,1,1,1.2)

postscript(file="crosssections.eps",width=8,height=10,horiz=F)
layout(t(matrix(1:12,nrow=3)),widt,hght)
plotno = 0
for(rs in condplot)
{
   plotno = plotno + 1
   if(((plotno-1) %% 3) == 0)
   {
      par(mar=c(0,4,0,0))
      par(xaxt="n")
      par(yaxt="s")
   } else
   {
      par(mar=c(0,0,0,0))
      par(xaxt="n")
      par(yaxt="n")
   }
   if(plotno > 8)
   {
      par(mar=c(4,0,0,0))
      par(xaxt="s")
      par(yaxt="n")
   }
   if(plotno == 10)

   {
      par(mar=c(4,4,0,0))
      par(xaxt="s")
      par(yaxt="s")
   }
   if(plotno == 11)
   {
      par(mar=c(4,0,0,0))
      par(xaxt="s")
      par(yaxt="n")
   }
   if(plotno == 9)
   {
      par(mar=c(0,0,0,0))
      par(xaxt="s")
      par(yaxt="n")
   }
   shft = log10(sum(datwghts))
   shft = shft - log10(dVc(rs)*1622/((180/pi)^2))
   plotconditional(results,2,varnames=c("Redshift","Absolute Magnitude"),val=rs,transtext="Logarithm of",trans=tran,normalize=F,
        reverse=T,lcol=1,ylimits=c(-10,-4),lwd=2,shft=shft,plotaxes=F,errcol="black")
   box()
   if(plotno > 8)
   {
      axis(1,at=c(-24,-25,-26,-27,-28,-29,-30),labels=c("","-25","","-27","","-29",""),cex.axis=1.5)
      mtext(expression(M),side=1,line=2.5,cex=1.1)
   }
   if(any(plotno == c(1,4,7,10)))
   {
      axis(2,at=c(-10,-9,-8,-7,-6,-5,-4),labels=c("","-9","","-7","","-5",""),cex.axis=1.5)
      mtext(expression(Phi(M)~~(Mpc^{-3}*mag^{-1})),side=2,line=2.5,cex=1.1)
   }
   text(-29.5,-4.5,paste("z =",as.character(rs)),cex=1.75)
   hold1 = RichLumfunc[RichLumfunc[,1]==rs,]
   pos = round(rs*100) + 1
   Richcol=rgb(0.65,0.65,0.65)
   points(hold1[,2],hold1[,3],lwd=2,lty=1,type="b",col=Richcol,pch=16)
   shftx = 0.1
   for(i in 1:length(hold1[,2]))
   {
      lo = log10(max(0.00000000000000001,(10^(hold1[i,3])+hold1[i,4]/(10^9))))
      hi = log10(max(0.00000000000000001,(10^(hold1[i,3])-hold1[i,4]/(10^9))))
      lines(rep(hold1[i,2],2),c(lo,hi),col=Richcol,lwd=2,lty=1)
      lines(rep(hold1[i,2],2)-c(-shftx,shftx),c(lo,lo),col=Richcol,lwd=2,lty=1)
      lines(rep(hold1[i,2],2)-c(-shftx,shftx),c(hi,hi),col=Richcol,lwd=2,lty=1)
   }
}
dev.off()