← Back to branch summary

~ubuntu-branches/ubuntu/gutsy/proj/gutsy

« back to all changes in this revision

Viewing changes to src/PJ_sconics.c

  • Committer: Bazaar Package Importer
  • Author(s): Peter S Galbraith
  • Date: 2002-01-11 10:27:12 UTC
  • Revision ID: james.westby@ubuntu.com-20020111102712-ayi18r8y2eesv0y9
Tags: upstream-4.4.5
Import upstream version 4.4.5

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/PJ_sconics.c
 
1
#ifndef lint
 
2
static const char SCCSID[]="@(#)PJ_sconics.c    4.1     94/05/22        GIE     REL";
 
3
#endif
 
4
#define PROJ_PARMS__ \
 
5
        double  n; \
 
6
        double  rho_c; \
 
7
        double  rho_0; \
 
8
        double  sig; \
 
9
        double  c1, c2; \
 
10
        int             type;
 
11
#define PJ_LIB__
 
12
#include        <projects.h>
 
13
#define EULER 0
 
14
#define MURD1 1
 
15
#define MURD2 2
 
16
#define MURD3 3
 
17
#define PCONIC 4
 
18
#define TISSOT 5
 
19
#define VITK1 6
 
20
#define EPS10   1.e-10
 
21
#define EPS 1e-10
 
22
#define LINE2 "\n\tConic, Sph\n\tlat_1= and lat_2="
 
23
PROJ_HEAD(tissot, "Tissot")
 
24
        LINE2;
 
25
PROJ_HEAD(murd1, "Murdoch I")
 
26
        LINE2;
 
27
PROJ_HEAD(murd2, "Murdoch II")
 
28
        LINE2;
 
29
PROJ_HEAD(murd3, "Murdoch III")
 
30
        LINE2;
 
31
PROJ_HEAD(euler, "Euler")
 
32
        LINE2;
 
33
PROJ_HEAD(pconic, "Perspective Conic")
 
34
        LINE2;
 
35
PROJ_HEAD(vitk1, "Vitkovsky I")
 
36
        LINE2;
 
37
/* get common factors for simple conics */
 
38
        static int
 
39
phi12(PJ *P, double *del) {
 
40
        double p1, p2;
 
41
        int err = 0;
 
42
 
 
43
        if (!pj_param(P->params, "tlat_1").i ||
 
44
                !pj_param(P->params, "tlat_2").i) {
 
45
                err = -41;
 
46
        } else {
 
47
                p1 = pj_param(P->params, "rlat_1").f;
 
48
                p2 = pj_param(P->params, "rlat_2").f;
 
49
                *del = 0.5 * (p2 - p1);
 
50
                P->sig = 0.5 * (p2 + p1);
 
51
                err = (fabs(*del) < EPS || fabs(P->sig) < EPS) ? -42 : 0;
 
52
                *del = *del;
 
53
        }
 
54
        return err;
 
55
}
 
56
FORWARD(s_forward); /* spheroid */
 
57
        double rho;
 
58
 
 
59
        switch (P->type) {
 
60
        case MURD2:
 
61
                rho = P->rho_c + tan(P->sig - lp.phi);
 
62
                break;
 
63
        case PCONIC:
 
64
                rho = P->c2 * (P->c1 - tan(lp.phi));
 
65
                break;
 
66
        default:
 
67
                rho = P->rho_c - lp.phi;
 
68
                break;
 
69
        }
 
70
        xy.x = rho * sin( lp.lam *= P->n );
 
71
        xy.y = P->rho_0 - rho * cos(lp.lam);
 
72
        return (xy);
 
73
}
 
74
INVERSE(s_inverse); /* ellipsoid & spheroid */
 
75
        double rho;
 
76
 
 
77
        rho = hypot(xy.x, xy.y = P->rho_0 - xy.y);
 
78
        if (P->n < 0.) {
 
79
                rho = - rho;
 
80
                xy.x = - xy.x;
 
81
                xy.y = - xy.y;
 
82
        }
 
83
        lp.lam = atan2(xy.x, xy.y) / P->n;
 
84
        switch (P->type) {
 
85
        case PCONIC:
 
86
                lp.phi = atan(P->c1 - rho / P->c2) + P->sig;
 
87
                break;
 
88
        case MURD2:
 
89
                lp.phi = P->sig - atan(rho - P->rho_c);
 
90
                break;
 
91
        default:
 
92
                lp.phi = P->rho_c - rho;
 
93
        }
 
94
        return (lp);
 
95
}
 
96
FREEUP; if (P) pj_dalloc(P); }
 
97
        static PJ *
 
98
setup(PJ *P) {
 
99
        double del, cs;
 
100
        int i;
 
101
 
 
102
        if( (i = phi12(P, &del)) )
 
103
                E_ERROR(i);
 
104
        switch (P->type) {
 
105
        case TISSOT:
 
106
                P->n = sin(P->sig);
 
107
                cs = cos(del);
 
108
                P->rho_c = P->n / cs + cs / P->n;
 
109
                P->rho_0 = sqrt((P->rho_c - 2 * sin(P->phi0))/P->n);
 
110
                break;
 
111
        case MURD1:
 
112
                P->rho_c = sin(del)/(del * tan(P->sig)) + P->sig;
 
113
                P->rho_0 = P->rho_c - P->phi0;
 
114
                P->n = sin(P->sig);
 
115
                break;
 
116
        case MURD2:
 
117
                P->rho_c = (cs = sqrt(cos(del))) / tan(P->sig);
 
118
                P->rho_0 = P->rho_c + tan(P->sig - P->phi0);
 
119
                P->n = sin(P->sig) * cs;
 
120
                break;
 
121
        case MURD3:
 
122
                P->rho_c = del / (tan(P->sig) * tan(del)) + P->sig;
 
123
                P->rho_0 = P->rho_c - P->phi0;
 
124
                P->n = sin(P->sig) * sin(del) * tan(del) / (del * del);
 
125
                break;
 
126
        case EULER:
 
127
                P->n = sin(P->sig) * sin(del) / del;
 
128
                del *= 0.5;
 
129
                P->rho_c = del / (tan(del) * tan(P->sig)) + P->sig;     
 
130
                P->rho_0 = P->rho_c - P->phi0;
 
131
                break;
 
132
        case PCONIC:
 
133
                P->n = sin(P->sig);
 
134
                P->c2 = cos(del);
 
135
                P->c1 = 1./tan(P->sig);
 
136
                if (fabs(del = P->phi0 - P->sig) - EPS10 >= HALFPI)
 
137
                        E_ERROR(-43);
 
138
                P->rho_0 = P->c2 * (P->c1 - tan(del));
 
139
                break;
 
140
        case VITK1:
 
141
                P->n = (cs = tan(del)) * sin(P->sig) / del;
 
142
                P->rho_c = del / (cs * tan(P->sig)) + P->sig;
 
143
                P->rho_0 = P->rho_c - P->phi0;
 
144
                break;
 
145
        }
 
146
        P->inv = s_inverse;
 
147
        P->fwd = s_forward;
 
148
        P->es = 0;
 
149
        return (P);
 
150
}
 
151
ENTRY0(euler) P->type = EULER; ENDENTRY(setup(P))
 
152
ENTRY0(tissot) P->type = TISSOT; ENDENTRY(setup(P))
 
153
ENTRY0(murd1) P->type = MURD1; ENDENTRY(setup(P))
 
154
ENTRY0(murd2) P->type = MURD2; ENDENTRY(setup(P))
 
155
ENTRY0(murd3) P->type = MURD3; ENDENTRY(setup(P))
 
156
ENTRY0(pconic) P->type = PCONIC; ENDENTRY(setup(P))
 
157
ENTRY0(vitk1) P->type = VITK1; ENDENTRY(setup(P))