# plot data, regression, lowess smooth
plot_data_and_fit = function(indata,fit,n,xunits,yunits,title,colfit,pch1,titlesize,textsize){
  # fit model
  #fit=lm(y~x+I(x^2)+sin(2*pi*x)+cos(2*pi*x),indata) #linear+quadratic+annual
  #fit=lm(y~x+sin(2*pi*x)+cos(2*pi*x),indata) #linear+annual
  #fit=lm(y~x,indata) #linear
  #slope=coef(summary(fit))[2,1]
  #slopeerror=2*coef(summary(fit))[2,2]#2 sigma
  fit=gls(y~x,data=indata,corr=corARMA(p=1,q=1))
  #fit=gls(y~x+sin(2*pi*x)+cos(2*pi*x),data=indata,corr=corARMA(p=1,q=1))
  #fit=gls(y~x+I(x^2),data=indata,corr=corARMA(p=1,q=1))
  fitsummmary=summary(fit)
  slope = fitsummmary$tTable[2,1]
  slopeerror = 2*fitsummmary$tTable[2,2]#2 sigma
  plot(indata$x,indata$y,type="o",pch=pch1,lwd=2,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=xunits,ylab=yunits,main=title)
  points(indata$x,fit$fit,type="l",lwd=2,lty=2,col=colfit)
  lowerbound=fit$fit-slopeerror*indata$x
  lowerbound=lowerbound - mean(lowerbound) + mean(fit$fit)
  points(indata$x,lowerbound,type="l",lwd=3,lty=1,col=colfit)
  upperbound=fit$fit+slopeerror*indata$x
  upperbound=upperbound - mean(upperbound) + mean(fit$fit)
  points(indata$x,upperbound,type="l",lwd=3,lty=1,col=colfit)
  confint(fit,digits=6)
  midpoint=(indata$x[n]+indata$x[1])/2.0
  top=(indata$y[which.max(indata$y)]-indata$y[which.min(indata$y)])*0.99+indata$y[which.min(indata$y)]
  text(midpoint,top,sprintf("%+.3f ± %.3f %s/%s",slope,slopeerror,yunits,xunits),cex=2,col=colfit)
  #par(new=T)
  #plot(lowess(indata$y,indata$x), col="blue",axes=F,xlab="",ylab="") # lowess smooth
}

# fit trends at different timespans
fit_many_trends = function(indata,n,stepsize,tlims,xunits,yunits,title,colfit,pch1,titlesize,textsize){
  #ntrend = n-30 #so the smallest subset has enough points
  ntrend = n-10 #so the smallest subset has enough points
  backtimes = indata$x[1:ntrend] #hold fits from different starting points to the end.
  backtrends = indata$y[1:ntrend] #hold trend fits from different starting points to the end.
  backtrenderrors = indata$y[1:ntrend] #hold 2 sigma error bars on trend fits from different starting points to the end.
  #Loop in reverse so smallest subset is first because it sometimes produces an error.
  for(j in seq(ntrend,1,by=-stepsize)){
    cdata = indata[j:n,]
    # remove mean again
    cdata$y = cdata$y - mean(cdata$y)
    #tfit=lm(y~x+sin(2*pi*x)+cos(2*pi*x),cdata)
    #tfit=lm(y~x,cdata)
    #backtrends[j] = coef(summary(tfit))[2,1]
    #backtrenderrors[j] = 2*coef(summary(tfit))[2,2]#2 sigma
    tfit=gls(y~x,data=cdata,corr=corARMA(p=1,q=1))
    tsummmary=summary(tfit)
    backtrends[j] = tsummmary$tTable[2,1]
    backtrenderrors[j] = 2*tsummmary$tTable[2,2]#2 sigma
  }
  if(stepsize>1){
    backtimes = backtimes[seq(ntrend,1,by=-stepsize)]
    backtrends = backtrends[seq(ntrend,1,by=-stepsize)]
    backtrenderrors = backtrenderrors[seq(ntrend,1,by=-stepsize)]
  }
  ttitle=paste(title,"trend to",sep=" ")
  ttitle=paste(ttitle,sprintf("%.1f",indata$x[n]),sep=" ")
  txunits=paste("starting",xunits,sep=" ")
  tyunits=paste(yunits,xunits,sep="/")
  if(length(tlims)==1) tlims=c(min(backtrends),max(backtrends))
  plot(backtimes,backtrends,type="o",pch=pch1,lwd=2,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=txunits,ylab=tyunits,main=ttitle,ylim=tlims)
  abline(h=0)
  upperbound=backtrends+backtrenderrors
  lowerbound=backtrends-backtrenderrors
  points(backtimes,upperbound,type="l",pch=pch1,lwd=1,col=colfit)
  points(backtimes,lowerbound,type="l",pch=pch1,lwd=1,col=colfit)
  cat('Trends end on ',indata$x[n],'.\n',sep='')
  for(j in length(backtimes):1) {
    cat('Trend starting on ',backtimes[j],' has trend ',backtrends[j],' +/- ',backtrenderrors[j],' ',yunits,'/',xunits,'.\n',sep='')
  }
}

# fit accels at different timespans
fit_many_accels = function(indata,n,stepsize,alims,xunits,yunits,title,colfit,pch1,titlesize,textsize){
  #naccel = n-90 #so the smallest subset has enough points
  naccel = n-10 #so the smallest subset has enough points
  backtimes = indata$x[1:naccel] #hold fits from different starting points to the end.
  backaccels = indata$y[1:naccel] #hold accel fits from different starting points to the end.
  backaccelerrors = indata$y[1:naccel] #hold 2 sigma error bars on accel fits from different starting points to the end.
  #Loop in reverse so smallest subset is first because it sometimes produces an error.
  for(j in seq(naccel,1,by=-stepsize)){
    cdata = indata[j:n,]
    # remove mean again
    cdata$y = cdata$y - mean(cdata$y)
    #afit=lm(y~x+I(x^2)+sin(2*pi*x)+cos(2*pi*x),cdata)
    #afit=lm(y~x+I(x^2),cdata)
    #Note that accel is twice the quadratic coefficient.
    #backaccels[j] = 2*coef(summary(afit))[3,1]
    #backaccelerrors[j] = 4*coef(summary(afit))[3,2]#2 sigma
    afit=gls(y~x+I(x^2),data=cdata,corr=corARMA(p=1,q=1))
    asummmary=summary(afit)
    #Note that accel is twice the quadratic coefficient.
    backaccels[j] = 2*asummmary$tTable[3,1]
    backaccelerrors[j] = 4*asummmary$tTable[3,2]#2 sigma
  }
  if(stepsize>1){
    backtimes = backtimes[seq(naccel,1,by=-stepsize)]
    backaccels = backaccels[seq(naccel,1,by=-stepsize)]
    backaccelerrors = backaccelerrors[seq(naccel,1,by=-stepsize)]
  }
  summary(afit)
  atitle=paste(title,"accel. to",sep=" ")
  atitle=paste(atitle,sprintf("%.1f",indata$x[n]),sep=" ")
  axunits=paste("starting",xunits,sep=" ")
  ayunits=paste(yunits,xunits,sep="/")
  ayunits=paste(ayunits,"^2")
  if(length(alims)==1) alims=c(min(backaccels),max(backaccels))
  plot(backtimes,backaccels,type="o",pch=pch1,lwd=2,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=axunits,ylab=ayunits,main=atitle,ylim=alims)
  abline(h=0)
  upperbound=backaccels+backaccelerrors
  lowerbound=backaccels-backaccelerrors
  points(backtimes,upperbound,type="l",pch=pch1,lwd=1,col=colfit)
  points(backtimes,lowerbound,type="l",pch=pch1,lwd=1,col=colfit)
}

# piecewise fit
piecewise_fit = function(indata,n,xunits,yunits,title,col1,col2,pch1,titlesize,textsize){
  library(segmented)
  linearfit=lm(y~x,indata)
  #linearfit=lm(y~x+sin(2*pi*x)+cos(2*pi*x),indata) #linear+annual
  midpoint=(indata$x[n]+indata$x[1])/2.0
  segments <- segmented(linearfit, seg.Z = ~x, psi=midpoint)
  summary(segments)
  slope(segments)
  #Why does the second std. error not match the value given in slope()? The first does!
  coef(summary(segments))
  slope1=coef(summary(segments))[2,1]
  slope2=coef(summary(segments))[3,1]
  slope2=slope1+slope2#second slope is sum of the two coefficients
  confint(segments,digits=6)
  plot(indata$x,indata$y,type="o",pch=pch1,lwd=2,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=xunits,ylab=yunits,main=title)
  #plot the segments and their slopes
  plot(segments,add=T,lwd=3,col=c(col1,col2))
  top=(indata$y[which.max(indata$y)]-indata$y[which.min(indata$y)])*0.99+indata$y[which.min(indata$y)]
  bottom=(indata$y[which.max(indata$y)]-indata$y[which.min(indata$y)])*0.03+indata$y[which.min(indata$y)]
  text(midpoint,top,sprintf("%.3f %s/%s",slope1,yunits,xunits),cex=2,col=col1)
  text(midpoint,bottom,sprintf("%.3f %s/%s",slope2,yunits,xunits),cex=2,col=col2)
  #plot the breakpoint and its 95% confidence interval
  lines(segments,col=3,pch=19,bottom=T,lwd=3)
  #draw.history(segments)
  #plot(segments,type="l",col="blue",add=T)
}

# plot de-annualized(?) data
plot_deannualized_data = function(indata,n,xunits,yunits,title,colfit,pch1,titlesize,textsize){
  seasons=lm(y~sin(2*pi*x)+cos(2*pi*x),indata)
  dmass=indata$y-seasons$fit
  dtitle=paste("De-annualized",title,sep=" ")
  plot(indata$x,indata$y,pch=pch1,type="o",lwd=2,cex.main=titlesize,cex.axis=textsize,cex.lab=textsize,xlab=xunits,ylab=yunits,main=dtitle)
  points(indata$x,dmass,col=colfit,pch=pch1,type="l") # regression
}