fork download
copy 
  1. #include <stdio.h>
  2. #include <math.h>
  3.  
  4. #define M_PI 3.141592653589793
  5.  
  6. void gsi_gauss_kruger(double lat, double lon, double lat0, double lon0, double k0, double *x, double *y) { 
  7. 	/* GRS80 楕円体定数 */
  8. 	const double a = 6378137.0;
  9. 	const double f = 1.0 / 298.257222101;
  10. 	const double e2 = 2.0 * f - f * f;
  11. 	const double e4 = e2 * e2;
  12. 	const double e6 = e4 * e2;
  13.  
  14. 	/* 子午線弧長 M(φ)（国土地理院の近似式） */
  15. 	const double A0 = 1.0 - e2/4.0 - 3.0*e4/64.0 - 5.0*e6/256.0;
  16. 	const double A2 = 3.0/8.0 * (e2 + e4/4.0 + 15.0*e6/128.0);
  17. 	const double A4 = 15.0/256.0 * (e4 + 3.0*e6/4.0);
  18. 	const double A6 = 35.0 * e6 / 3072.0;
  19.  
  20. 	/* 度 → ラジアン */ 
  21. 	double phi = lat * M_PI / 180;
  22. 	double phi0 = lat0 * M_PI / 180;
  23.  
  24. 	double sinp = sin(phi);
  25. 	double cosp = cos(phi);
  26. 	double tanp = tan(phi);
  27.  
  28. 	double M = a * (A0*phi - A2*sin(2*phi) + A4*sin(4*phi) - A6*sin(6*phi));
  29. 	double M0 = a * (A0*phi0 - A2*sin(2*phi0) + A4*sin(4*phi0) - A6*sin(6*phi0));
  30. 	double N = a / sqrt(1.0 - e2 * sinp * sinp);	/* 卯酉線曲率半径 */
  31. 	double eta2 = e2 * cosp * cosp / (1.0 - e2);	/* η^2 = e'^2 cos^2 φ */
  32.  
  33. 	/* 近似式の級数展開 */ 
  34. 	double A = cosp * (lon - lon0) * M_PI / 180;	/* 経度差 */
  35. 	double A2p = A*A;
  36. 	double A3p = A2p*A;
  37. 	double A4p = A2p*A2p;
  38. 	double A5p = A4p*A;
  39. 	double A6p = A4p*A2p;
  40.  
  41. 	/* 北向き x */ 
  42. 	*x = k0 * (M - M0 + N * tanp * ( A2p/2.0 + (5.0 - tanp*tanp + 9.0*eta2 + 4.0*eta2*eta2) * A4p / 24.0 + (61.0 - 58.0*tanp*tanp + tanp*tanp*tanp*tanp + 270.0*eta2 - 330.0*tanp*tanp*eta2) * A6p / 720.0 ));
  43.  
  44. 	/* 東向き y */
  45. 	*y = k0 * N * ( A + (1.0 - tanp*tanp + eta2) * A3p / 6.0 + (5.0 - 18.0*tanp*tanp + tanp*tanp*tanp*tanp + 14.0*eta2 - 58.0*tanp*tanp*eta2) * A5p / 120.0 );
  46. }
  47.  
  48. int main(void) {
  49. 	double x, y;
  50.  
  51. 	gsi_gauss_kruger(36.103774791666666, 140.08785504166664, 36.0, 139.8333333333333, 0.9999, &x, &y);
  52.  
  53. 	printf("x=%lf, y=%lf\n", x, y);
  54. 	return 0;
  55. }
  56.  
Success #stdin #stdout 0s 5320KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
x=11543.688323, y=22916.243554
https://ideone.com/cixKhk
language:
C (gcc 8.3)
created:
66 days ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!