Changeset 8698


Ignore:
Timestamp:
06/23/11 11:01:20 (14 years ago)
Author:
Mathieu Morlighem
Message:

some cosmetics

Location:
issm/trunk/src/c/modules
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/modules/Ll2xyx/Ll2xyx.cpp

    r8689 r8698  
    33
    44#include "./Ll2xyx.h"
     5#include "../../include/include.h"
    56#include "../../shared/shared.h"
    6 #include "../../include/include.h"
    7 #include "../../toolkits/toolkits.h"
    8 #include "../../EnumDefinitions/EnumDefinitions.h"
     7#include <math.h>
    98
    10 int Ll2xyx(double* x, double* y,
    11                    double* lat, double* lon, int ncoord,
    12                    int sgn){
    13 
     9int Ll2xyx(double* x, double* y, double* lat, double* lon, int ncoord, int sgn){
    1410/*  This is a cpp conversion of the following:
    1511%LL2XY - converts lat long to polar stereographic
     
    2622%                               -1 : south latitude (default is mer=0  lat=71)
    2723*/
    28 
    2924        double  delta,slat;
    3025        bool    flag=true;
    3126
    3227        /*  Get central_meridian and standard_parallel depending on hemisphere  */
    33         if      (sgn ==  1) {
    34                 delta = 45;
     28        if (sgn ==  1) {
     29                delta= 45;
    3530                slat = 70;
    3631                _printf_(flag,"Info: creating coordinates in polar stereographic (Std Latitude: 70N Meridian: 45).\n");
    3732        }
    3833        else if (sgn == -1) {
    39                 delta = 0;
     34                delta= 0;
    4035                slat = 71;
    4136                _printf_(flag,"Info: creating coordinates in polar stereographic (Std Latitude: 71S Meridian: 0).\n");
    4237        }
    43         else
    44                 _error_("Sign should be either +1 or -1.\n");
     38        else _error_("Sign should be either +1 or -1.\n");
    4539
    46         return(Ll2xyx(x,y,
    47                                   lat,lon,ncoord,
    48                                   sgn,delta,slat));
     40        return(Ll2xyx(x,y, lat,lon,ncoord, sgn,delta,slat));
    4941}
    5042
    51 int Ll2xyx(double* x, double* y,
    52                    double* lat, double* lon, int ncoord,
    53                    int sgn, double central_meridian, double standard_parallel){
    54 
     43int Ll2xyx(double* x, double* y, double* lat, double* lon, int ncoord, int sgn, double central_meridian, double standard_parallel){
    5544/*  This is a cpp conversion of the following:
    5645%LL2XY - converts lat long to polar stereographic
     
    7463        double  T,rho,sl,tc,mc;
    7564
     65<<<<<<< .mine
     66        if((sgn!=1) && (sgn!=-1)) _error_("Sign should be either +1 or -1");
     67=======
    7668        if      ((sgn !=  1) && (sgn != -1))
    7769                _error_("Sign should be either +1 or -1.\n");
     70>>>>>>> .r8690
    7871
    7972        delta = central_meridian;
     
    8174
    8275        /*  Conversion constant from degrees to radians  */
    83         cde  = 57.29577951;
     76        cde = 57.29577951;
    8477        /*  Radius of the earth in meters  */
    85         re   = 6378.273*pow(10.,3.);
     78        re  = 6378.273*pow(10.,3.);
    8679        /*  Eccentricity of the Hughes ellipsoid squared  */
    87         ex2   = .006693883;
     80        ex2 = 0.006693883;
    8881        /*  Eccentricity of the Hughes ellipsoid  */
    89         ex    =  sqrt(ex2);
    90         /*  pi is not defined in cpp  */
    91         pi   = 4.*atan(1.);
     82        ex  =  sqrt(ex2);
    9283
    9384        /*  loop over all the coordinate pairs  */
    94 
    95         for (i=0; i<ncoord; i++) {
    96                 latitude  = fabs(lat[i]) * pi/180.;
    97                 longitude = (lon[i] + delta) * pi/180.;
     85        for(i=0; i<ncoord; i++){
     86                latitude  = fabs(lat[i]) * PI/180.;
     87                longitude = (lon[i] + delta) * PI/180.;
    9888
    9989                /*  compute X and Y in grid coordinates.  */
    100                 T = tan(pi/4.-latitude/2.) / pow(((1.-ex*sin(latitude))/(1.+ex*sin(latitude))),(ex/2.));
     90                T = tan(PI/4.-latitude/2.) / pow(((1.-ex*sin(latitude))/(1.+ex*sin(latitude))),(ex/2.));
    10191
    102                 if ((90. - slat) <  1.e-5)
     92                if ((90. - slat) < 1.e-5)
    10393                        rho = 2.*re*T/sqrt(pow((1.+ex),(1.+ex))*pow((1.-ex),(1.-ex)));
    10494                else {
    105                         sl  = slat*pi/180.;
    106                         tc  = tan(pi/4.-sl/2.)/pow(((1.-ex*sin(sl))/(1.+ex*sin(sl))),(ex/2.));
     95                        sl  = slat*PI/180.;
     96                        tc  = tan(PI/4.-sl/2.)/pow(((1.-ex*sin(sl))/(1.+ex*sin(sl))),(ex/2.));
    10797                        mc  = cos(sl)/sqrt(1.0-ex2*(pow(sin(sl),2.)));
    10898                        rho = re*mc*T/tc;
    10999                }
    110100
    111                 y[i] = -rho * (double)sgn * cos((double)sgn*longitude);
    112                 x[i] =  rho * (double)sgn * sin((double)sgn*longitude);
     101                y[i]= -rho*(double)sgn*cos(sgn*longitude);
     102                x[i]=  rho*(double)sgn*sin(sgn*longitude);
    113103
    114                 if (latitude >= pi / 2.) {
     104                if (latitude>= PI/2.){
    115105                        x[i] = 0.0;
    116106                        y[i] = 0.0;
     
    118108                }
    119109        }
    120 
    121110        return(iret);
    122111}
    123 
  • issm/trunk/src/c/modules/Ll2xyx/Ll2xyx.h

    r8673 r8698  
    77
    88/* local prototypes: */
    9 int Ll2xyx(double* x, double* y,
    10                    double* lat, double* lon, int ncoord,
    11                    int sgn);
    12 
    13 int Ll2xyx(double* x, double* y,
    14                    double* lat, double* lon, int ncoord,
    15                    int sgn, double central_meridian, double standard_parallel);
     9int Ll2xyx(double* x, double* y, double* lat, double* lon, int ncoord, int sgn);
     10int Ll2xyx(double* x, double* y, double* lat, double* lon, int ncoord, int sgn, double central_meridian, double standard_parallel);
    1611
    1712#endif  /* _LL2XYX_H */
    18 
  • issm/trunk/src/c/modules/Xy2llx/Xy2llx.cpp

    r8689 r8698  
    33
    44#include "./Xy2llx.h"
     5#include "../../include/include.h"
    56#include "../../shared/shared.h"
    6 #include "../../include/include.h"
    7 #include "../../toolkits/toolkits.h"
    8 #include "../../EnumDefinitions/EnumDefinitions.h"
     7#include <math.h>
    98
    10 int Xy2llx(double* lat, double* lon,
    11                    double* x, double* y, int ncoord,
    12                    int sgn){
    13 
     9int Xy2llx(double* lat, double* lon, double* x, double* y, int ncoord, int sgn){
    1410/*  This is a cpp conversion of the following:
    1511%XY2LL - converts xy to lat long
     
    3228
    3329        /*  Get central_meridian and standard_parallel depending on hemisphere  */
    34         if      (sgn == 1) {
    35                 delta = 45;
     30        if (sgn == 1) {
     31                delta= 45;
    3632                slat = 70;
    3733                _printf_(flag,"Warning: expecting coordinates in polar stereographic (Std Latitude: 70N Meridian: 45).\n");
    3834        }
    3935        else if (sgn == -1) {
    40                 delta = 0;
     36                delta= 0;
    4137                slat = 71;
    4238                _printf_(flag,"Warning: expecting coordinates in polar stereographic (Std Latitude: 71S Meridian: 0).\n");
    4339        }
    44         else
    45                 _error_("Sign should be either +1 or -1.\n");
     40        else _error_("Sign should be either +1 or -1.\n");
    4641
    47         return(Xy2llx(lat,lon,
    48                                   x,y,ncoord,
    49                                   sgn,delta,slat));
     42        return(Xy2llx(lat,lon,x,y,ncoord,sgn,delta,slat));
    5043}
    5144
    52 int Xy2llx(double* lat, double* lon,
    53                    double* x, double* y, int ncoord,
    54                    int sgn, double central_meridian, double standard_parallel){
    55 
     45int Xy2llx(double* lat, double* lon, double* x, double* y, int ncoord, int sgn, double central_meridian, double standard_parallel){
    5646/*  This is a cpp conversion of the following:
    5747%XY2LL - converts xy to lat long
     
    7565        double  sl,rho,cm,T,chi;
    7666
    77         if      ((sgn !=  1) && (sgn != -1))
    78                 _error_("Sign should be either +1 or -1.\n");
     67        if((sgn!=1) && (sgn!=-1)) _error_("Sign should be either +1 or -1.\n");
    7968
    8069        delta = central_meridian;
     
    8675        re   = 6378.273*pow(10.,3.);
    8776        /*  Eccentricity of the Hughes ellipsoid squared  */
    88         ex2   = .006693883;
     77        ex2  = 0.006693883;
    8978        /*  Eccentricity of the Hughes ellipsoid  */
    90         ex    =  sqrt(ex2);
    91         /*  pi is not defined in cpp  */
    92         pi   = 4.*atan(1.);
     79        ex   =  sqrt(ex2);
    9380
    9481        /*  loop over all the coordinate pairs  */
     82        for(i=0; i<ncoord; i++){
     83                sl = slat*PI/180.;
     84                cm = cos(sl)/sqrt(1.0-ex2*(pow(sin(sl),2.)));
     85                rho= sqrt(pow(x[i],2.) + pow(y[i],2.));
     86                T  = tan((PI/4.0) - (sl/2.0))/pow(((1.0-ex*sin(sl))/(1.0+ex*sin(sl))),(ex/2.0));
    9587
    96         for (i=0; i<ncoord; i++) {
    97                 sl  = slat*pi/180.;
    98                 cm = cos(sl) / sqrt(1.0 - ex2 * (pow(sin(sl),2.)));
    99                 rho = sqrt(pow(x[i],2.) + pow(y[i],2.));
    100                 T = tan((pi / 4.0) - (sl / 2.0)) / pow(((1.0 - ex * sin(sl)) / (1.0 + ex * sin(sl))),(ex / 2.0));
     88                if(fabs(slat-90.) < 1.e-5)
     89                        T =rho*sqrt(pow((1.+ex),(1.+ex))*pow((1.-ex),(1.-ex)))/2./re;
     90                else
     91                        T =rho*T/(re*cm);
    10192
    102                 if  (fabs(slat-90.) < 1.e-5)
    103                         T = rho * sqrt(pow((1. + ex),(1. + ex)) * pow((1. - ex),(1. - ex))) / 2. / re;
    104                 else
    105                         T = rho * T / (re * cm);
    106 
    107                 chi = (pi / 2.0) - 2.0 * atan(T);
     93                chi = (PI / 2.0) - 2.0 * atan(T);
    10894                lat[i] = chi + ((ex2 / 2.0) + (5.0 * pow(ex2,2.0) / 24.0) + (pow(ex2,3.0) / 12.0)) *
    10995                           sin(2.0 * chi) + ((7.0 * pow(ex2,2.0) / 48.0) + (29.0 * pow(ex2,3.0) / 240.0)) *
     
    114100                lon[i] = (double)sgn * lon[i];
    115101
    116                 if (rho <= 0.1) {
     102                if(rho <= 0.1){
    117103                        lat[i] = 90. * (double)sgn;
    118104                        lon[i] = 0.0;
     
    127113        return(iret);
    128114}
    129 
  • issm/trunk/src/c/modules/Xy2llx/Xy2llx.h

    r8673 r8698  
    77
    88/* local prototypes: */
    9 int Xy2llx(double* lat, double* lon,
    10                    double* x, double* y, int ncoord,
    11                    int sgn);
    12 
    13 int Xy2llx(double* lat, double* lon,
    14                    double* x, double* y, int ncoord,
    15                    int sgn, double central_meridian, double standard_parallel);
     9int Xy2llx(double* lat, double* lon, double* x, double* y, int ncoord, int sgn);
     10int Xy2llx(double* lat, double* lon, double* x, double* y, int ncoord, int sgn, double central_meridian, double standard_parallel);
    1611
    1712#endif  /* _XY2LLX_H */
    18 
Note: See TracChangeset for help on using the changeset viewer.