###CALCULATE TWO NEXT NEAREST STATIONS AND INSERT IN DETAILS TABLE	

##################
###  FUNCTIONS
##################
		
### 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
     }		

###Make rural network for given id 	
   pi180=pi/180;  
   nextnearest_mann=function(i) {
	   temp_rur=!is.na(info$mann)
	   station_versions=info[temp_rur,] ;dim(station_versions)
	   station_versions$dist=NA
	   lat0=details$lat[i];long0=details$long[i]
	   temp1= (abs (station_versions$lat-lat0) <20) #
		#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(station_versions$long -long0)
	   delta= apply( cbind(x %%360,(360-x)%%360),1,min) 
  	   temp2=(delta<test_angle) 	#second screen: within 1000 km by longitude
	   if(sum(temp1&temp2)<2) temp2=rep(TRUE,1732)
	   ruralstations=station_versions[temp1&temp2,]
		#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
	   order1=order(ruralstations$dist,ruralstations$mannid)
	   nextnearest_mann=	 ruralstations$jones[order1][1:2]
	   nextnearest_mann
	 }


	download.file("http://www.climateaudit.org/data/mann.2008/details.tab","temp.dat",mode="Wb")
	load("temp.dat")
	details$next1=NA;details$next2=NA
	for (i in 1:1209) details[i,c("next1","next2")]=nextnearest_mann(i)