| 1 | // %pacpus:license{
|
|---|
| 2 | // This file is part of the PACPUS framework distributed under the
|
|---|
| 3 | // CECILL-C License, Version 1.0.
|
|---|
| 4 | // %pacpus:license}
|
|---|
| 5 | /// @file
|
|---|
| 6 | /// @author Firstname Surname <firstname.surname@utc.fr>
|
|---|
| 7 | /// @date Month, Year
|
|---|
| 8 | /// @version $Id: transf_lamb_93.h 76 2013-01-10 17:05:10Z kurdejma $
|
|---|
| 9 | /// @copyright Copyright (c) UTC/CNRS Heudiasyc 2006 - 2013. All rights reserved.
|
|---|
| 10 | /// @brief Brief description.
|
|---|
| 11 | ///
|
|---|
| 12 | /// Detailed description.
|
|---|
| 13 |
|
|---|
| 14 | #include <math.h>
|
|---|
| 15 | #define pi 3.141592653589793
|
|---|
| 16 | double atanh(double z)
|
|---|
| 17 | {
|
|---|
| 18 | return (0.5 * log((z+1)/(1-z)));
|
|---|
| 19 | }
|
|---|
| 20 |
|
|---|
| 21 | double abs_me(double x)
|
|---|
| 22 | {
|
|---|
| 23 | if(x<0) x=-x;
|
|---|
| 24 | return x;
|
|---|
| 25 | }
|
|---|
| 26 |
|
|---|
| 27 | //transformation des longitudes tatitudes en Lambert93
|
|---|
| 28 | void lonlattolam(double lon, double lat, double & lam93x, double & lam93y)
|
|---|
| 29 | {
|
|---|
| 30 | double GRS_a = 6378137;
|
|---|
| 31 | double GRS_f = 1./298.257222101;
|
|---|
| 32 |
|
|---|
| 33 | double GRS_b = GRS_a*(1-GRS_f);
|
|---|
| 34 | double GRS_bb= GRS_b*GRS_b;
|
|---|
| 35 | double GRS_aa= 40680631590769.;
|
|---|
| 36 | double GRS_e = sqrt((GRS_aa - GRS_bb) / (GRS_aa));
|
|---|
| 37 |
|
|---|
| 38 | double n = 0.725607765053267;
|
|---|
| 39 | double C = 11754255.4261;
|
|---|
| 40 | double XS = 700000;
|
|---|
| 41 | double YS = 12655612.0499;
|
|---|
| 42 |
|
|---|
| 43 | double latiso;
|
|---|
| 44 | latiso = atanh(sin(lat)) - GRS_e*atanh(GRS_e*sin(lat));
|
|---|
| 45 | double gamma;
|
|---|
| 46 | gamma = (lon - 0.0523598775598299)*n;
|
|---|
| 47 | double R;
|
|---|
| 48 | R = C * exp(-n*latiso);
|
|---|
| 49 |
|
|---|
| 50 | lam93x = R *sin(gamma)+XS;
|
|---|
| 51 | lam93y = -R *cos(gamma)+YS;
|
|---|
| 52 | //printf("x: %.20lf\n y: %.20lf\n",lam93x,lam93y);
|
|---|
| 53 |
|
|---|
| 54 | }
|
|---|
| 55 |
|
|---|
| 56 |
|
|---|
| 57 |
|
|---|
| 58 | //Transformation des coordonnees de Lambert93 en Longitude lattitude
|
|---|
| 59 |
|
|---|
| 60 | void lamtolonlat(double lamx, double lamy, double & lon, double & lat)
|
|---|
| 61 | {
|
|---|
| 62 |
|
|---|
| 63 | double GRS_a = 6378137;
|
|---|
| 64 | double GRS_f = 1./298.257222101;
|
|---|
| 65 | double GRS_b = GRS_a*(1-GRS_f);
|
|---|
| 66 | double GRS_bb= GRS_b*GRS_b;
|
|---|
| 67 | double GRS_aa= 40680631590769.;
|
|---|
| 68 | double GRS_e = sqrt((GRS_aa - GRS_bb) / (GRS_aa));
|
|---|
| 69 |
|
|---|
| 70 | //double n = 0.725607765053267;
|
|---|
| 71 | //double C = 11754255.4261;
|
|---|
| 72 | //double XS = 700000;
|
|---|
| 73 | //double YS = 12655612.0499;
|
|---|
| 74 |
|
|---|
| 75 |
|
|---|
| 76 | //lamx = lamx-700000;
|
|---|
| 77 | //lamy = lamy-12655612.0499;
|
|---|
| 78 |
|
|---|
| 79 | double gamma;
|
|---|
| 80 | gamma = atan(-(lamx-700000)/(lamy-12655612.0499));
|
|---|
| 81 |
|
|---|
| 82 |
|
|---|
| 83 | lon = gamma/0.725607765053267 + 0.0523598775598299;
|
|---|
| 84 | double R;
|
|---|
| 85 | R = sqrt((lamx-700000) * (lamx-700000) + (lamy-12655612.0499) * (lamy-12655612.0499));
|
|---|
| 86 |
|
|---|
| 87 | double latiso;
|
|---|
| 88 | latiso = log((11754255.4261)/R)/(0.725607765053267);
|
|---|
| 89 |
|
|---|
| 90 | double phiNew, phiOld;
|
|---|
| 91 | phiOld =1;
|
|---|
| 92 | phiNew= asin (tanh ( latiso + GRS_e * atanh(GRS_e * sin(phiOld))));
|
|---|
| 93 | //printf("\nphiNew: %.20lf",phiNew);
|
|---|
| 94 | while (abs_me(phiOld-phiNew) > 1e-10)
|
|---|
| 95 | {
|
|---|
| 96 |
|
|---|
| 97 | if(abs_me(phiOld-phiNew) > 1e-10)
|
|---|
| 98 | {
|
|---|
| 99 |
|
|---|
| 100 | phiOld = phiNew;
|
|---|
| 101 | phiNew = asin(tanh(latiso+GRS_e*atanh(GRS_e*sin(phiOld))));
|
|---|
| 102 | }
|
|---|
| 103 | else
|
|---|
| 104 | phiOld = phiNew;
|
|---|
| 105 | }
|
|---|
| 106 |
|
|---|
| 107 | lat = phiNew;
|
|---|
| 108 |
|
|---|
| 109 | }
|
|---|
