Initial commit

1 files changed, 200 insertions, 0 deletions

diff --git a/vo/votools/gasplib/ccgsxy.x b/vo/votools/gasplib/ccgsxy.x
new file mode 100644
index 00000000..588dacef
--- /dev/null
+++ b/vo/votools/gasplib/ccgsxy.x
@@ -0,0 +1,200 @@
+define	MAX_ITERATIONS		20
+# CCGSXY --  Routine for computing X,Y pixel position from a given RA and Dec
+#            on a GSSS plate.
+
+procedure ccgsxy (plate_centre_ra,
+     		  plate_centre_dec,
+     		  plate_centre_x,
+     		  plate_centre_y,
+     		  x_pixel_size,
+     		  y_pixel_size,
+     		  plate_scale,
+     		  amd_x,
+     		  amd_y,              
+     		  object_x,
+     		  object_y,
+     		  object_mag,
+     		  object_col,
+     		  object_ra,
+     		  object_dec)
+
+
+double	plate_centre_ra		#Plate Right Ascension (radians)
+double	plate_centre_dec	#Plate Declination     (radians)
+double	plate_centre_x		#Position used in soln.(microns)
+double	plate_centre_y		# "
+double	x_pixel_size		#Scan pixel size (microns)
+double	y_pixel_size		# "
+double	plate_scale		#Plate scale (arcsec/mm)
+double	amd_x[20]		#Plate model coefficients
+double	amd_y[20]		# "
+double	object_x		#Pixel position for object
+double	object_y		# "
+double	object_mag		#Object magnitude 
+double	object_col		#Colour
+double	object_ra		#Object RA (radians)	
+double	object_dec		# "     Dec
+double	xi_object		#Standard coords (arcsec)
+double	eta_object		# "
+
+double	x			#Position from centre (mm)
+double	y			# "
+double	delta_x			#correction to x
+double	delta_y			#correction to y
+int	n_iterations		#no.iterations used
+
+double	f			#Function for x model
+double	g			#Function for y model
+double	fx			#Deriv. x model wrt x
+double	gx			#Deriv. y model wrt x
+double	fy			#Deriv. x model wrt y
+double	gy			#deriv. y model wrt y
+
+int	ierr			#Stats flag
+
+begin
+
+	#
+	# Convert RA and Dec to standard coordinates
+	#
+	call treqst (plate_centre_ra, plate_centre_dec,
+     	             object_ra, object_dec, xi_object, eta_object)
+	
+	# Iterate by Newtons method
+	n_iterations = 0
+	ierr = 0
+	
+	# Get initial estimate of x,y
+	x = xi_object / plate_scale
+	y = eta_object / plate_scale
+
+    while (ierr == 0) {
+	n_iterations=n_iterations+1
+
+	#
+	# X plate model
+	#
+	f = amd_x(1)*x 			     	   	+
+     	    amd_x(2)*y 		   	   	   	+
+     	    amd_x(3)   		       	    	   	+
+     	    amd_x(4)*x*x 	       	    	   	+ 
+     	    amd_x(5)*x*y	       	      	   	+
+     	    amd_x(6)*y*y	       	      	   	+
+     	    amd_x(7)*(x*x+y*y)	       		   	+
+     	    amd_x(8)*x*x*x	       		   	+
+     	    amd_x(9)*x*x*y	       	       	   	+
+     	    amd_x(10)*x*y*y	       	       	   	+
+     	    amd_x(11)*y*y*y	       	       	   	+
+     	    amd_x(12)*x*(x*x+y*y)   		   	+
+     	    amd_x(13)*x*(x*x+y*y)**2	       	   	+
+     	    amd_x(14)*object_mag		   	+
+     	    amd_x(15)*object_mag**2		   	+
+     	    amd_x(16)*object_mag**3		   	+
+     	    amd_x(17)*object_mag*x		   	+
+     	    amd_x(18)*object_mag*(x*x+y*y)		+
+      	    amd_x(19)*object_mag*x*(x*x+y*y)	       	+ 
+     	    amd_x(20)*object_col			-
+     	              xi_object
+	#
+	# Derivative of X model wrt x
+	#
+	fx = amd_x(1)	 			       	+
+     		amd_x(4)*2.0*x 		    		       	+ 
+     		amd_x(5)*y	       			       	+       	
+     		amd_x(7)*2.0*x				       	+
+     		amd_x(8)*3.0*x*x     			       	+
+     		amd_x(9)*2.0*x*y			       	+
+     		amd_x(10)*y*y		       		       	+
+     		amd_x(12)*(3.0*x*x+y*y)	      		       	+
+     		amd_x(13)*(5.0*x**4+6.0*x*x*y*y+y**4)  	       	+
+     		amd_x(17)*object_mag		       	       	+
+     		amd_x(18)*object_mag*2.0*x	       	       	+
+      		amd_x(19)*object_mag*(3.0*x*x+y*y)	
+	#
+	# Derivative of X model wrt y
+	#
+	fy = amd_x(2)	 		      	       	+
+     		amd_x(5)*x		       		       	+
+     		amd_x(6)*2.0*y		       		       	+
+     		amd_x(7)*2.0*y		  		       	+
+     		amd_x(9)*x*x	   			       	+
+     		amd_x(10)*x*2.0*y			       	+
+     		amd_x(11)*3.0*y*y  			       	+
+     		amd_x(12)*2.0*x*y		      	       	+
+     		amd_x(13)*4.0*x*y*(x*x+y*y)		       	+
+     		amd_x(18)*object_mag*2.0*y		       	+
+      		amd_x(19)*object_mag*2.0*x*y		
+	#
+	# Y plate model
+	#
+	g =	amd_y(1)*y  				       	+
+     		amd_y(2)*x 		       		       	+
+     		amd_y(3)   		   		       	+
+     		amd_y(4)*y*y 	     			       	+ 
+     		amd_y(5)*y*x	      			       	+
+     		amd_y(6)*x*x	      			       	+
+     		amd_y(7)*(x*x+y*y)     			       	+
+     		amd_y(8)*y*y*y				       	+
+     		amd_y(9)*y*y*x		   		       	+
+     		amd_y(10)*y*x*x		 		       	+
+     		amd_y(11)*x*x*x		  		       	+
+     		amd_y(12)*y*(x*x+y*y)	  		       	+
+     		amd_y(13)*y*(x*x+y*y)**2	    	       	+
+     		amd_y(14)*object_mag			       	+
+     		amd_y(15)*object_mag**2			       	+
+     		amd_y(16)*object_mag**3		  	       	+
+     		amd_y(17)*object_mag*y		  	       	+
+     		amd_y(18)*object_mag*(x*x+y*y)	   	       	+
+      		amd_y(19)*object_mag*y*(x*x+y*y)	       	+
+     		amd_y(20)*object_col		      	       	-
+     		          eta_object
+	#
+	# Derivative of Y model wrt x
+	#
+	gx = amd_y(2)	 			       	+
+     		amd_y(5)*y		 		       	+
+     		amd_y(6)*2.0*x	    			       	+
+     		amd_y(7)*2.0*x	      			       	+
+     		amd_y(9)*y*y				       	+
+     		amd_y(10)*y*2.0*x			       	+
+     		amd_y(11)*3.0*x*x		   	       	+
+     		amd_y(12)*2.0*x*y		  	       	+
+     		amd_y(13)*4.0*x*y*(x*x+y*y)  		       	+
+     		amd_y(18)*object_mag*2.0*x		       	+
+      		amd_y(19)*object_mag*y*2.0*x		
+	#
+	# Derivative of Y model wrt y
+	#
+	gy = amd_y(1)				       	+
+     		amd_y(4)*2.0*y 		   	       	       	+ 
+     		amd_y(5)*x				       	+
+     		amd_y(7)*2.0*y		     		       	+
+     		amd_y(8)*3.0*y*y	      		       	+
+     		amd_y(9)*2.0*y*x		       	       	+
+     		amd_y(10)*x*x		     		       	+
+     		amd_y(12)*3.0*y*y		       	       	+
+     		amd_y(13)*(5.0*y**4+6.0*x*x*y*y+x**4) 	       	+
+     		amd_y(17)*object_mag			       	+
+     		amd_y(18)*object_mag*2.0*y	  	       	+
+      		amd_y(19)*object_mag*(x*x+3.0*y*y)		
+
+	delta_x = (-f * gy + g * fy) / (fx * gy - fy * gx)
+	delta_y = (-g * fx + f * gx) / (fx * gy - fy * gx)
+	x = x + delta_x
+	y = y + delta_y
+
+	if (dmax1 (dabs(delta_x), dabs(delta_y),
+     	           dabs(f), dabs(g)) < 1.e-5) 
+           ierr=1
+
+	if (n_iterations == MAX_ITERATIONS) 
+	   ierr=2
+
+     } 
+	 
+	# Convert x,y from mm about plate centre to pixels
+
+	object_x = (plate_centre_x - x*1000.0d0) / x_pixel_size
+	object_y = (plate_centre_y + y*1000.0d0) / y_pixel_size
+
+end