options(echo = FALSE);
options(digits=15);
####################################################################
## Execute following on command prompt (omit the "##"):
## /usr/bin/R -q --no-restore --no-save < dzp_vs_tempfit.r
##
####################################################################

library(MASS);

# returns the nearest occurence of x in vec
nearest.vec <- function(x, vec) {
  smallCandidate <- findInterval(x, sort(vec), all.inside=TRUE);
  largeCandidate <- smallCandidate + 1;
  #nudge is TRUE if large candidate is nearer, FALSE otherwise
  nudge <- 2 * x > vec[smallCandidate] + vec[largeCandidate];
  return(smallCandidate + nudge);
}

#----------------------------------------------------------------
# Input filename (from Sherry)
# cols: scan   temp    w1dmag   w1sig   w2dmag   w2sig

infile1="popsv7_scan_temp_dmags_hdr.txt";
outfile1="dzp_vs_temp_int.v7.5.pdf";
outfile2="dzp_vs_temp_int.v7.5_residuals.pdf";
#----------------------------------------------------------------

  data<-scan(infile1,list("",0,0,0,0,0));
  scan<-data[[1]];
  mjd<-data[[2]];
  w1dmag<-data[[3]];
  w1sig<-data[[4]];
  w2dmag<-data[[5]];
  w2sig<-data[[6]];
  
  pdf(file=outfile1, pointsize=12);
  par(oma=c(2,2,2,2));
  plot(mjd, w1dmag, type="p", pch=20, cex=0.7,
       xlab="",
       ylab="w?dmag [mag]",
       ylim=c(-0.02,0.01),
       col="blue",cex.lab=1.3,cex.axis=1.3,cex.main=1.3,
       xaxt="n");
  points(mjd, w2dmag, col="green", pch=20, cex=0.7);
  abline(h=0,col="black",lty=2);

  # X-axis (bottom).
  mjdlab <- seq(min(mjd),max(mjd),by=(max(mjd)-min(mjd))/20);
  mjdstr <- sprintf("%5.2f",mjdlab);
  axis(1, las=2, at=mjdlab, labels=mjdstr);
  mtext(1, text="bsplmout temperature [K]", line=5, cex=1.3);

  # X-axis (top).
  idx<-nearest.vec(mjdlab,mjd);
  scanlab <- scan[idx];
  axis(3, las=2, at=mjdlab, labels=scanlab);
  mtext(3, text="approx. scan id", line=4.3, cex=1.3);

  # robust linear fits.
  lin <- rlm(w1dmag ~ mjd + 1, method="MM");
  par <- coefficients(summary(lin));
  intercept <- par[1,1];
  slope <- par[2,1];
  intercept; slope;
  abline(intercept, slope, col="blue");
  w1resid <- w1dmag - (slope*mjd + intercept);

  i<-which((mjd < 73.465) | (mjd > 73.56));
  lin <- rlm(w2dmag[i] ~ mjd[i] + 1, method="MM");
  par <- coefficients(summary(lin));
  intercept <- par[1,1];
  slope <- par[2,1];
  intercept; slope;
  abline(intercept, slope, col="green");
  w2resid <- w2dmag - (slope*mjd + intercept);

  legend(mjdlab[[2]], 0.01, c("W1",
                              "W2"),
         col = c("blue",
                 "green"),
         text.col="black",
         pch = c(20,20),
         merge=FALSE, bty='n');

#----------------------------------------------------------------
# plot residuals in above fits.

  pdf(file=outfile2, pointsize=12);
  plot(mjd, w1resid, type="p", pch=4, cex=1.2,
       xlab="bsplmout temperature [K]",
       ylab="obs - predicted w?dmag",
       ylim=c(-0.01,0.01),
       col="blue",cex.lab=1.3,cex.axis=1.3,cex.main=1.3);
  points(mjd, w2resid, col="green", pch=4, cex=1.2);
  abline(h=0,col="black",lty=2);

  sd(w1resid); sd(w2resid);
  median(w1resid); median(w2resid[i]);

q();