Stata code to transform GPS to Rijksdriehoekscoordinates

/****************************************/
/* GPS TO RIJKSDRIEHOEKSCOORDINATEN		*/
/****************************************/

// IMPORTANT: Make sure that the variables "rijks_x" & "rijks_y" don't exist yes
// Drop "rijks_x" & "rijks_y"
*capture drop rijks_x
*capture drop rijks_y

/****************************************/
/* PART 1: CALCULATION MODULE			*/
/****************************************/

// Program has two arguments (parameters) as input
// Command syntax: <program name> <latitude to be transformed> <name of longitude variable>
// Example: gps2rijks_calc_module 53.160753042 4.824761912
capture program drop gps2rijks_calc_module	// Check if program exists and if it does, drop it
program gps2rijks_calc_module
quietly {
	// Constants
	local x0 = 155000
	local y0 = 463000
	local phi0 = 52.1551729
	local lam0 = 5.387203657

	// GPS input -> change to variable input
	local cord_x = `1'		// test value = 53.160753042 should result in rijks_x = 117380.11
	local cord_y = `2'		// test value = 4.824761912 should result in rijks_y = 575040.4

	// Delta variables
	local dPhi = 0.36*(`cord_x' - `phi0')
	local dLam = 0.36*(`cord_y' - `lam0')

	// Fill arrays with calculation parameters
	mata :	Rp = 0,1,2,0,1,3,1,0,2
	mata :	Rq = 1,1,1,3,0,1,3,2,3
	mata :	Rpq = 190094.945,-11832.228,-114.221,-32.391,-0.705,-2.340,-0.608,-0.008,0.148
	mata :	Sp = 1,0,2,1,3,0,2,1,0,1
	mata :	Sq = 0,2,0,2,0,1,2,1,4,4
	mata :	Spq = 309056.544,3638.893,73.077,-157.984,59.788,0.433,-6.439,-0.032,0.092,-0.054

	// Initiate/reset variables for calculation loops
	local x = 0
	local y = 0
	local i = 1
	local j = 1

	// Calculate X Rijksdriekhoekscoordinaat
	while `i' <= 9 {
		mata : st_local("x", strofreal(`x' + (Rpq[1,`i']) * (`dPhi')^Rp[1,`i']* (`dLam')^Rq[1,`i']))
		local ++i
	}
	global res_rijks_x = `x' + `x0'	// Save rijksdriehoekscoordinaat x to global

	// Calculate Y Rijksdriekhoekscoordinaat
	while `j' <= 10 {
		mata : st_local("y", strofreal(`y' + (Spq[1,`j']) * (`dPhi')^Sp[1,`j']* (`dLam')^Sq[1,`j']))
		local ++j
	}
	global res_rijks_y = `y' + `y0'	// Save rijksdriehoekscoordinaat y to global
}
end


/****************************************/
/* PART 2: INPUT & SAVE MODULE			*/
/****************************************/

// Program has two arguments (parameters) as input
// Command syntax: <program name> <name of latitude variable> <name of longitude variable>
// Example: gps2rijks latitude longitude
capture program drop gps2rijks	// Check if program exists and if it does, drop it
program gps2rijks
quietly {
	// Capture latitude and longitude variable names
	local lat = "`1'"
	local lon = "`2'"

	// Generate empty variables
	gen rijks_x = .
	gen rijks_y = .

	// Replace empty variable values with the calculated Rijksdriehoekscoordinaten
	local i = 1
	while `i' <= _N {
		gps2rijks_calc_module `lat'[`i'] `lon'[`i'] // Call calculation module with GPS coordinates of the i'th observation
		replace rijks_x = $res_rijks_x if _n == `i'
		replace rijks_y = $res_rijks_y if _n == `i'
		local `++i'
	}
}
end