# Functions for Tripole grid
#  Murray, R. J., 1996: Explicit Generation of Orthogonal Grids for Ocean Models.
#  J. Comp. Phys., 126, 251-273.

RADEG = 180. / pi
DTOR  = pi / 180.

# This function computes grid cell edges that will 
# have edges at the poles and an even number around the
# globe.

Choose_tripole_longitudes <- function(
		tripole_longitude,
		dlongitude,
		inDegrees ){
		
		dlon = dlongitude
		if ( inDegrees ) dlon = dlon * DTOR
		
		intervals_per_pi = round( pi / dlon ) 
		intervals = 2. * intervals_per_pi
		dlon = pi / intervals_per_pi
		if ( inDegrees ) dlon = dlon * RADEG
		
		x=tripole_longitude + dlon * c(0:intervals)
		
}

# The great circle distance between two points, returned in Radians
# This algorithm is better if the arc is small.

Great_circle_distance <- function(lon1,lat1,lon2,lat2){

    2.0*asin(sqrt((sin((lat1-lat2)*0.5))^2 + 
                 cos(lat1)*cos(lat2)*(sin((lon1-lon2)*0.5))^2)) 
}



Global_Tripole <- function( 
       longitudes,			# Longitudes to use in the spherical region
       latitudes,			# vector of deisred latitudes along tripole symmetry meridian
       tripole_longitude,	#   longitude of "first" tripole 
       tripole_latitude,	#   latitude of polar cap tripoles
       inDegrees=FALSE,		#	if TRUE, input is in degrees
       name=NULL			#   if present, use as an ID
       )
{

	gsw = G_tripole( longitudes, latitudes,
			tripole_longitude, tripole_latitude, inDegrees=inDegrees )
	gc  = Find_grid_centers( gsw )
	gw  = Find_W_faces( gsw )
	gs  = Find_S_faces( gsw )
	
	im = gsw$im
	jm = gsw$jm
	j3 = gsw$jcap1
	j2 = j3-1
	
	# we now know the grid box corners and centers
	# Now compute the grid face lengths. These are
	# approximate, in that they presume that a great 
	# circle distance can be used to determine the length
		
	dx_S = matrix( nrow=im-1, ncol=jm )
	dy_W = matrix( nrow=im,   ncol=jm-1 )
		
		dx_S[1:(im-1),1:j2] = cos( gsw$y[1:(im-1),1:j2] )*
							(gsw$x[2:im,1:j2]-gsw$x[1:(im-1),1:j2])
		
		# Since we know the corners PLUS the centers of each grid box,
		# we can compute the distance better by using 2 great circle
		# arcs for each side.  The first arc goes from the corner to the
		# midpoint, while the second goes from the mid-point to the other
		# corner.
		
		dx_S[1:(im-1),j3:jm] = Great_circle_distance( 
					lon1=gsw$x[1:(im-1),j3:jm],
					lat1=gsw$y[1:(im-1),j3:jm],
					lon2=gs$x[1:(im-1),j3:jm],
					lat2=gs$y[1:(im-1),j3:jm] ) + 
				Great_circle_distance(
					lon1=gs$x[1:(im-1),j3:jm],
					lat1=gs$y[1:(im-1),j3:jm],
					lon2=gsw$x[2:im,j3:jm],
					lat2=gsw$y[2:im,j3:jm] )  
				
				
		dy_W[1:im,1:(j2-1)] = gsw$y[1:im,2:j2]-gsw$y[1:im,1:(j2-1)]
		dy_W[1:im,j2:(jm-1)] = Great_circle_distance( 
					lon1=gsw$x[1:im,j2:(jm-1)],
					lat1=gsw$y[1:im,j2:(jm-1)],
					lon2=gw$x[1:im,j2:(jm-1)],
					lat2=gw$y[1:im,j2:(jm-1)] ) + 
				Great_circle_distance( 
					lon1=gw$x[1:im,j2:(jm-1)],
					lat1=gw$y[1:im,j2:(jm-1)],
					lon2=gsw$x[1:im,j3:jm],
					lat2=gsw$y[1:im,j3:jm] )
	
		dx_c = matrix( nrow=im-1, ncol = jm-1 )
		dy_c = matrix( nrow=im-1, ncol = jm-1 )
		
		# The distance between the grid faces is dx_c
		# we will also compute it as two arcs:  one from
		# the face to the grid center, the other from the
		# grid center to the opposite face
		
		dx_c[1:(im-1),1:j2] = cos( gc$y[1:(im-1),1:j2] )*
							(gw$x[2:im,1:j2]-gw$x[1:(im-1),1:j2])
		
		# Since we know the corners PLUS the centers of each grid box,
		# we can compute the distance better by using 2 great circle
		# arcs for each side.  The first arc goes from the corner to the
		# midpoint, while the second goes from the mid-point to the other
		# corner.
		
		dx_c[1:(im-1),j3:(jm-1)] = Great_circle_distance( 
					lon1=gw$x[1:(im-1),j3:(jm-1)],
					lat1=gw$y[1:(im-1),j3:(jm-1)],
					lon2=gc$x[1:(im-1),j3:(jm-1)],
					lat2=gc$y[1:(im-1),j3:(jm-1)] ) + 
				Great_circle_distance(
					lon1=gc$x[1:(im-1),j3:(jm-1)],
					lat1=gc$y[1:(im-1),j3:(jm-1)],
					lon2=gw$x[2:im,j3:(jm-1)],
					lat2=gw$y[2:im,j3:(jm-1)] )  
				
				
		dy_c[1:(im-1),1:(j2-1)] = gs$y[1:(im-1),2:j2]-gs$y[1:(im-1),1:(j2-1)]
		dy_c[1:(im-1),j2:(jm-1)] = Great_circle_distance( 
					lon1=gs$x[1:(im-1),j2:(jm-1)],
					lat1=gs$y[1:(im-1),j2:(jm-1)],
					lon2=gc$x[1:(im-1),j2:(jm-1)],
					lat2=gc$y[1:(im-1),j2:(jm-1)] ) + 
				Great_circle_distance( 
					lon1=gc$x[1:(im-1),j2:(jm-1)],
					lat1=gc$y[1:(im-1),j2:(jm-1)],
					lon2=gs$x[1:(im-1),j3:jm],
					lat2=gs$y[1:(im-1),j3:jm] )
		
		Area = matrix(nrow=im-1, ncol=jm-1)
		
				Area[1:(im-1),1:(j2-1)] =
					dy_c[1:(im-1),1:(j2-1)] *
					dx_c[1:(im-1),1:(j2-1)]
				
				Area[1:(im-1),j2:(jm-1)] =
					(dy_W[1:(im-1),j2:(jm-1)] + dy_W[2:im,j2:(jm-1)])*
					(dx_S[1:(im-1),j2:(jm-1)] + dx_S[1:(im-1),(j2+1):jm])*
					0.25

		SWcorners = list( x=gsw$x, y=gsw$y, chi=gsw$chi )
		Centers   = list( x=gc$x , y=gc$y, chi=gc$chi, dx=dx_c, dy=dy_c, area=Area   )
		Sfaces 	  = list( x=gs$x , y=gs$y, chi=gs$chi, dx=dx_S   )
		Wfaces    = list( x=gw$x , y=gw$y, chi=gw$chi, dy=dy_W   )
		Inputs    = list( longitudes=gsw$longitudes, latitudes=gsw$latitudes,
		                  tripole_longitude=gsw$tripole_longitude,
		                  tripole_latitude=gsw$tripole_latitude,
		                  )
		                  
		                  
		if ( is.null(name) )  {
			Name = paste("Tripole_",(im-1),"x",(jm-1),"_(",
					tripole_longitude,",",tripole_latitude,")",
			    sep="")
		} else { Name=name }
		
		
		rv = list( im=im-1, jm=jm-1, SWcorners=SWcorners, Centers=Centers, 
				   Sfaces=Sfaces, Wfaces=Wfaces,
				   Inputs=Inputs,
				   Name = Name)
				   
		
}

G_tripole <- function( 
       longitudes,			# Longitudes to use in the spherical region
       latitudes,			# vector of deisred latitudes along tripole symmetry meridian
       tripole_longitude,	#   longitude of "first" tripole 
       tripole_latitude,	#   latitude of polar cap tripoles
       inDegrees=FALSE		#	if TRUE, input and output are in degrees
       )
{
	
	# Tripole_latitude  - latitude where the Tripole region should be matched
	# Tripole_longitude - longitude of the tripole closest to i=1 grid point

	
	if ( inDegrees ) {
		lons = longitudes * DTOR
		lats = latitudes  * DTOR
		Phi_P = tripole_latitude * DTOR
		Lam_P = tripole_longitude * DTOR
	} else {
		lons = longitudes
		lats = latitudes 
		Phi_P = tripole_latitude
		Lam_P = tripole_longitude
	}

	# Test some of the inputs:
	
	if ( lons[1] < Lam_P ) stop("First longitude must be >= tripole_longitude")


# Declare some functions here 

#  Earth_latitude - returns the longitude on the sphere
#  of a tri-polar grid with the tangent plane at the North pole

Earth_longitude <- function(Lambda_P,Grid_longitude,Grid_latitude){

		Lambda_P - atan2(  sin(Grid_longitude) , tan(Grid_latitude) )
}

Earth_latitude <- function(Grid_longitude,Grid_latitude){
	
		Xc= acos( cos(Grid_longitude) * cos(Grid_latitude) )
		pi*0.5 - 2.*atan( rP * tan(Xc*0.5) )
}

Choose_Xc <- function( desired_lats, Phi_P ) 2. * atan(  tan( (pi*0.5-desired_lats)*0.5 ) / rP ) 


Rotation_angle <- function( Lambda_P,lon,lat){

    
	ye1=Earth_latitude( lon,lat)
	
	coscos = cos(lon) * cos(lat)
	tanlat = tan(lat)
	
	{ dxdY = - cos(ye1) * cos(lon) * tanlat /
	            ( sin(lon)^2 + tanlat^2  ) }
	
	tanacos2 = tan( acos(coscos) * 0.5 )^2
	
	{ dydY =  - rP * ( 1. + tanacos2 ) * sin(lon) * cos(lat) /
	   ( ( sqrt( 1. - coscos^2 ) ) * ( 1. + rP^2 * tanacos2 ) ) }
	
    
    atan2(dxdY,-dydY)
    
}


	# Start by computing the actual latitude of the tripole
	#  We do not check consistency that the pole latitude is a multiple 
	# of the grid, because we might want the pole at a grid cell center point 
	# while the grid is computed for the grid cell corners.  But be careful!
	
	
	j1 = 1
	j2 = max( which( lats < Phi_P ) )
	j3 = j2 + 1
	jm = length(lats)
	
	
	# Now compute the diameter ratio of the grid sphere to the Earth sphere
	
	rP = tan( (pi*0.5 - Phi_P)*0.5 )

	
	# how many points in the polar cap?
	  
	cap_lats = lats[ j3:jm ]
	ncaps = length(cap_lats) 
	
	# Transform the desired latitudes along the symmetry meridian
	# into grid sphere longitudes
	
	lon_g = Choose_Xc( cap_lats, Phi_P )
	  
	# Where are points around cap?
	
		# check on longitudes[1] and tripole_longitude		
		  
	  	i2 = max( which( ( lons - Lam_P ) <= pi ))
	  
		xg1 = matrix( lon_g , i2,   ncaps, byrow=TRUE )
		yg1 = matrix( -pi*0.5 + lons[1:i2] - Lam_P , i2,   ncaps)
		ye1_cap=Earth_latitude( xg1,yg1)
		xe1_cap=Earth_longitude(Lam_P+pi,xg1,yg1)
		chi1_cap=Rotation_angle( Lam_P+pi,xg1,yg1 ) 
	  
		# Now do the second quadrant:	
		
		i3 = i2 + 1
		im = length(lons)
		n2 = im - i3 + 1
	  
		xg2 = matrix( -lon_g , n2, ncaps, byrow=TRUE )
		yg2 = matrix(  pi*0.5 - ( lons[i3:im] - Lam_P - pi ),  n2, ncaps)
		ye2_cap=Earth_latitude( xg2,yg2)
		xe2_cap=Earth_longitude(  Lam_P,xg2,yg2 ) + pi
		chi2_cap =Rotation_angle( Lam_P,xg2,yg2 ) + pi

		# Now patch the two quadrants into a single matrix
		#  rbind does this for us.
	  
		xe_cap  = rbind( xe1_cap, xe2_cap )
		ye_cap  = rbind( ye1_cap, ye2_cap )
		chi_cap = rbind( chi1_cap, chi2_cap )

		xe_south  = matrix( lons, im, j2 )
		ye_south  = matrix( lats[1:j2], im, j2, byrow=TRUE)
		chi_south = matrix( 0., im,j2)
	  
	  	# Now patch the sphere and cap together
		x   = cbind( xe_south, xe_cap )
		y   = cbind( ye_south, ye_cap )
		chi = cbind( chi_south, chi_cap )
		
		chi[ which( chi > pi ) ] = chi[ which( chi > pi ) ] - 2 * pi
		

		list( x=x, y=y, chi=chi, name="tripolar global grid",
	      longitudes=lons,latitudes=lats,
	      tripole_longitude=Lam_P,tripole_latitude=Phi_P,
	      inDegrees=inDegrees,im=im,jm=jm,jcap1=j3)
	
}

Find_grid_centers <- function( g ){

	lon_c = ( g$longitudes[1:(g$im-1)] + g$longitudes[2:g$im] ) * 0.5
	lat_c = ( g$latitudes[1:(g$jm-1)]  + g$latitudes[2:g$jm]  ) * 0.5
	
	G_tripole(lon_c,lat_c,
			g$tripole_longitude,g$tripole_latitude)
	
}

Find_W_faces <- function( g ){

	lat_c = ( g$latitudes[1:(g$jm-1)]  + g$latitudes[2:g$jm]  ) * 0.5
	
	G_tripole(g$longitudes[1:g$im],lat_c,
			g$tripole_longitude,g$tripole_latitude)
	
}

Find_S_faces <- function( g ){

	lon_c = ( g$longitudes[1:(g$im-1)] + g$longitudes[2:g$im] ) * 0.5
	
	G_tripole(lon_c,g$latitudes[1:g$jm],
			g$tripole_longitude,g$tripole_latitude)
	
}

WriteTripole <- function( g, filename ){

	con=file(filename,"wb")
	
	#	rv = list( im=im-1, jm=jm-1, SWcorners=SWcorners, Centers=Centers, 
	#			   Sfaces=Sfaces, Wfaces=Wfaces,
	#			   Inputs=Inputs,
	#			   Name = Name)
				   
	
	# First write out the SW corners
	writeBin(  c(g$im+1,g$jm+1,3), con )
	writeBin( as.vector( g$SWcorners$x ),  con )
	writeBin( as.vector( g$SWcorners$y ),  con )
	writeBin( as.vector( g$SWcorners$chi), con )
	
	writeBin(  c(g$im,g$j,6), con )
	writeBin( as.vector( g$Centers$x ),    con )
	writeBin( as.vector( g$Centers$y ),    con )
	writeBin( as.vector( g$Centers$chi ),  con )
	writeBin( as.vector( g$Centers$dx ),   con )
	writeBin( as.vector( g$Centers$dy ),   con )
	writeBin( as.vector( g$Centers$area ), con )
	
	writeBin( c(g$im,g$jm+1,4), con )
	writeBin( g$Sfaces$x,   con )
	writeBin( g$Sfaces$y,   con )
	writeBin( g$Sfaces$chi, con )
	writeBin( g$Sfaces$dx,  con )
	
	writeBin( c(g$im+1,g$jm,4), con )
	writeBin( g$Wfaces$x,   con )
	writeBin( g$Wfaces$y,   con )
	writeBin( g$Wfaces$chi, con )
	writeBin( g$Wfaces$dx,  con )
	
	writeBin( c(0,0,0), con )

	# 
	close(con)
}

plotTripole <- function( g , xskip=4, yskip=2, fill_color="seagreen" ) {

	im = g$im
	jm = g$jm

	j1 = min( which(g$SWcorners$y[1,1:jm+1] > 0 ) )

	map('world',interior=FALSE,proj="orthographic",fill=TRUE,col=fill_color)
	
	for (i in seq(1,im+1,xskip) ){ 
		mp=mapproject(g$SWcorners$x[i,j1:(jm+1)]*RADEG,g$SWcorners$y[i,j1:(jm+1)]*RADEG)
		lines( mp$x,mp$y ) 
	}
	for ( j in seq(j1,(jm+1),yskip) ){
		mp=mapproject(g$SWcorners$x[1:(im+1),j]*RADEG,g$SWcorners$y[1:(im+1),j]*RADEG)
		lines( mp$x,mp$y ) 
	}
}


tripole_example <- function(){

	lons = Choose_tripole_longitudes( -110., 10., inDegrees=TRUE)

	lats = seq( -90, 90, by=5. )
	
	g=Global_Tripole(lons,lats,-110.,60.,inDegrees=TRUE,name="10x5 Tripole")
	
	plotTripole(g,xskip=1,yskip=1); title('Tripole Grid',sub='GMU Poseidon Ocean Model: tripole.r \t\t\t             ref: Murray(1996)')
	
	WriteTripole( g, filename='tripole.dat' )
	
}