trim=function (chron) {  
	M0=range(time(chron)[!is.na(chron)])
	trim=ts(chron[ (time(chron)>=M0[1])&(time(chron)<=M0[2]) ],start=M0[1])
	trim}

long=function(x){ x=trim(x); M0=range(time(x));long=M0[2]-M0[1]+1;long}

hansen_delta=function(x,reference){
	combine=ts.union(x,reference)
	temp=!is.na(combine[,1])&!is.na(combine[,2])
	m0=apply(combine[temp,],2,mean)
	hansen_delta=m0[1]-m0[2]
	hansen_delta}

ssq=function(x){ssq=sum(x^2,na.rm=TRUE);ssq}

twoleg=function(x) {
	N=length(x)
	stat=rep(NA,N)
	temp=!is.na(x)
	index=(1:N)[temp]
	for (i in  index[2:(length(index)-1)]) {
	  z=data.frame(x, abs( (1:N)-i), (1:N)-i)   ;names(z)=c("y","abs","x")	
	  fm=lm(y~abs+x,data=z)
	  stat[i]=ssq(fm$residuals)	
	  }
	index0=(1:N)[(stat==min(stat,na.rm=TRUE))&!is.na(stat)]
     z=data.frame(x, abs( (1:N)-index0), (1:N)-i)   ;names(z)=c("y","abs","x")	
     twoleg=lm(y~abs+x,data=z)
     twoleg}

   ### circledist
	#calculates great circle distance  with earth radius 6378.137 km
	#from http://en.wikipedia.org/wiki/Great-circle_distance 2nd formula	
	
   circledist =function(loc, lat,long,R=6372.795) {
	N=length(lat)	
	if(N==0) circledist=NA else {
	  pi180=pi/180;  
   	  x= abs(long -loc[2])
	  if (N>1) delta= apply( cbind(x %%360,(360-x)%%360),1,min)  *pi180  else delta= min (x %%360,(360-x)%%360) *pi180  
	  loc=loc*pi180; lat=lat*pi180; long=long*pi180
	  theta= 2* asin( sqrt(  sin( (lat- loc[1])/2 )^2 + cos(lat)*cos(loc[1])* (sin(delta/2))^2       ))
	  circledist=R*theta
	  }
	circledist
     }		
	
	ruralstations=function(id0) {
	  if(is.numeric(id0)) i=id0 else i= (1:7364)[!is.na(match (names(giss.dset1),id0))] 
	   temp_rur=(stations$urban=="R")
	   stations$dist=NA
	   lat0=stations$lat[i];long0=stations$long[i]
	   temp1= (abs (stations$lat-lat0) <10) 
		#initial screening only by latitude: must be within 10 deg latitude to be within 1000 km
	   test_angle= 10/sin(abs(lat0)*pi180)
	   x=	abs(stations$long -long0)
	   delta= apply( cbind(x %%360,(360-x)%%360),1,min) 
  	   temp2=(delta<test_angle) 
		#second screen: within 1000 km by longitude
	   ruralstations=stations[temp1&temp2&temp_rur,c("id","name","long","lat","altitude","start_raw","end_raw","start_adj","end_adj","dist","pop","urban","lights")]
		#matrix of stations meeting first two screen tests
   	   ruralstations$dist=circledist(loc=c(lat0,long0),lat=ruralstations$lat,long=ruralstations$long)
		#calculate circle distance for all stations in first screen
	   temp=(ruralstations$dist<=1000)&!is.na(ruralstations$dist)
	   ruralstations=ruralstations[temp,]	
		#restrict to stations within 1000 km
	   ruralstations=ruralstations[order(ruralstations[,"dist"]),]
		#order by increasing distance
	   ruralstations$weights= 1- ruralstations$dist/1000
	   ruralstations
	}


##FUNCTION TO EMULATE STEP 2 COMBINING
   #requires giss.dset1; stat_dist   
 
   emulate_dset2=function(id0,method="base") {
	if (is.numeric(id0) &(id0<8000)) i=id0 else  i=(1:7364)[(names(giss.dset1)==id0)]
	rural=ruralstations(i);	K=nrow(rural)
		#this gives list of rural station ids
   	weights0= rural$weights		# weights according to Hansen method
	ruralid=rural$id
	info=stations[i,]

  if(method=="base"){
	chron=NULL
	for (k in 1:K){
	  n=(1:7364)[names(giss.dset1)==ruralid[k]]
	  test1=giss.dset1[[n]][[1]]
	  fred=hansenq(test1,method="hansen") 
	  chron=ts.union(chron,trim( round(fred[,"ANN"],2)) )
	  }
 	dimnames(chron)[[2]]=rural$name
	          #this collects annual network of rural stations using Hansen interpolation

	dset1.monthly=ts( c(t(giss.dset1[[i]][[1]][,2:13])) ,start=c(giss.dset1[[i]][[1]][1,1],1),freq=12)
		#monthly ts in this method  
		#dset1=round(trim(hansenq(dset1.monthly,method="ts")[,5]),2)
	dset1= ts(giss.dset1[[i]][[1]][,18], start=giss.dset1[[i]][[1]][1,1])

	dset2.monthly=ts( c(t(giss.dset2[[i]][[1]][,2:13])) ,start=c(giss.dset2[[i]][[1]][1,1],1),freq=12)
	dset2= ts(giss.dset2[[i]][[1]][,18], start=giss.dset2[[i]][[1]][1,1])
      }

 if(method=="read"){
	chron=read.table("http://www.climateaudit.org/data/giss/step2/compare.dat",sep="\t",header=TRUE)
	chron=ts(chron[,2:ncol(chron)],start=chron[1,1])
	          #this collects annual network of rural stations using Hansen interpolation

	target=read.table("http://www.climateaudit.org/data/giss/step2/target.dat",sep="\t",header=TRUE)	
	target=ts(target[,2:3],start=target[1,1])
	dset1= target[,1]
	dset2= target[,2]
	dset1.monthly=NA
	dset2.monthly=NA
 	rural=read.table("http://www.climateaudit.org/data/giss/step2/rural.dat",sep="\t",header=TRUE)	
	weights0=rural$weights
     }

	count=ts(apply(!is.na(chron),1,sum),start=tsp(chron)[1])
	count.adj=count;count.adj[count<3]=NA
		#if less than 3 stations, not calculated
	M0=range(time(count)[!is.na(count.adj)]) 
		#range of available

	Y=ts.union(dset1,count.adj)
	temp1=(time(Y)>=M0[1])&(time(Y)<=M0[2])&!is.na(Y[,1])&!is.na(Y[,2])
	N3=sum( temp1);  
		#number of years with both dset1 has values and count>=3 in range with count>=3
	N3/(M0[2]-M0[1]+1) # fraction of range 

  #do iteration  
     #set up arrays	
     	N=nrow(chron)
	W=t(array(rep(weights0,N),dim=c(length(weights0),N) ) )
	W[is.na(chron)]=NA
		#matrix of weights: at least 3 stations and station available
	long0=apply(chron,2,long)
		#calculates length of the candidate rural series
      	order1=order(long0,decreasing=TRUE);index=NULL
		#long0[order1] # orders by decreasing length
	delta=rep(NA,K);
   	fixed=NULL; available= 1:K; 
	m0=mean(dset1,na.rm=T)	

     #start with the longest series
       	j0=order1[1]; 
	reference=chron[,j0] #picks longest series to insert
	delta[j0]=hansen_delta(chron[,j0], dset1);
		#calculates Hansen delta uses Step 1 methods
	reference=reference-delta[j0];
	reference.seq=reference
	weights1=rep(weights0[j0],length(reference));weights1[is.na(reference)]=NA
	N=nrow(chron);fixed=j0

     #sequentially adds in the rural series 
	for (k in 2:K) {
	  j0=order1[k] #chooses the series to be added
   	  delta[j0]=hansen_delta(chron[,j0], reference);#calculates Hansen delta
	  weights2=W[,j0]
	  weights1=apply(as.matrix(W[,fixed]),1,sum,na.rm=T)
     		#g=function(x) g= weighted.mean(x,,1),na.rm=TRUE) 
	  
	  Y=cbind(chron[,j0]-delta[j0],reference,reference); X=cbind(weights2,weights1)
	  Y[,3]=NA
	  for(j in 1:N) Y[j,3]=  weighted.mean(Y[j,1:2],X[j,] ,na.rm=T) 
	  fixed=c(fixed,j0)
	  reference=Y[,3]
	  reference.seq=cbind(reference.seq,reference)
	}
	M0=range( time(reference)[count>=3],na.rm=T)

	Y=ts.union(dset1,dset2,reference,reference,reference,count)
	dimnames(Y)[[2]]=c("dset1","dset2","reference","adjustment","adjusted","count")
	Y[,4:5]=NA
	temp1= (Y[,"count"]>=3)&!is.na(Y[,1])&!is.na(Y[,"count"]);sum(temp1)  #patch assumption
	M0=range( c(time(Y)) [temp1],na.rm=T)
	temp=(time(Y)>=M0[1])# &(time(Y)<=M0[2])  #patch assumption #by obs only start seems to matter
	#Y[!temp,3]=NA
	adj0=twoleg(Y[temp,3]-Y[temp,1])
	Y[temp&!is.na(Y[,1]),4]=round(fitted(adj0),2)
	Y[,5]=round(Y[,1]+Y[,4],2)

	emulation=Y

	emulate_dset2=list(dset1.monthly,dset1,dset2.monthly,dset2,chron,count,W,delta,reference,info,adj0,emulation)
	names(emulate_dset2)=c("dset1.monthly","dset1","dset2.monthly","dset2","chron","count","weights","delta","reference","info","twoleg","emulation") 
	emulate_dset2
       }