Bienvenue sur PostGIS.fr

Bienvenue sur PostGIS.fr , le site de la communauté des utilisateurs francophones de PostGIS.

PostGIS ajoute le support d'objets géographique à la base de données PostgreSQL. En effet, PostGIS "spatialise" le serverur PostgreSQL, ce qui permet de l'utiliser comme une base de données SIG.

Maintenu à jour, en fonction de nos disponibilités et des diverses sorties des outils que nous testons, nous vous proposons l'ensemble de nos travaux publiés en langue française.

source: trunk/workshop-routing-foss4g/web/proj4js/lib/proj4js-combined.js @ 76

Revision 76, 175.0 KB checked in by djay, 12 years ago (diff)

Ajout du répertoire web

  • Property svn:executable set to *
Line 
1/*
2  proj4js.js -- Javascript reprojection library.
3 
4  Authors:      Mike Adair madairATdmsolutions.ca
5                Richard Greenwood richATgreenwoodmap.com
6                Didier Richard didier.richardATign.fr
7                Stephen Irons
8  License:      LGPL as per: http://www.gnu.org/copyleft/lesser.html
9                Note: This program is an almost direct port of the C library
10                Proj4.
11*/
12/* ======================================================================
13    proj4js.js
14   ====================================================================== */
15
16/*
17Author:       Mike Adair madairATdmsolutions.ca
18              Richard Greenwood rich@greenwoodmap.com
19License:      LGPL as per: http://www.gnu.org/copyleft/lesser.html
20
21$Id: Proj.js 2956 2007-07-09 12:17:52Z steven $
22*/
23
24/**
25 * Namespace: Proj4js
26 *
27 * Proj4js is a JavaScript library to transform point coordinates from one
28 * coordinate system to another, including datum transformations.
29 *
30 * This library is a port of both the Proj.4 and GCTCP C libraries to JavaScript.
31 * Enabling these transformations in the browser allows geographic data stored
32 * in different projections to be combined in browser-based web mapping
33 * applications.
34 *
35 * Proj4js must have access to coordinate system initialization strings (which
36 * are the same as for PROJ.4 command line).  Thes can be included in your
37 * application using a <script> tag or Proj4js can load CS initialization
38 * strings from a local directory or a web service such as spatialreference.org.
39 *
40 * Similarly, Proj4js must have access to projection transform code.  These can
41 * be included individually using a <script> tag in your page, built into a
42 * custom build of Proj4js or loaded dynamically at run-time.  Using the
43 * -combined and -compressed versions of Proj4js includes all projection class
44 * code by default.
45 *
46 * Note that dynamic loading of defs and code happens ascynchrously, check the
47 * Proj.readyToUse flag before using the Proj object.  If the defs and code
48 * required by your application are loaded through script tags, dynamic loading
49 * is not required and the Proj object will be readyToUse on return from the
50 * constructor.
51 *
52 * All coordinates are handled as points which have a .x and a .y property
53 * which will be modified in place.
54 *
55 * Override Proj4js.reportError for output of alerts and warnings.
56 *
57 * See http://trac.osgeo.org/proj4js/wiki/UserGuide for full details.
58*/
59
60/**
61 * Global namespace object for Proj4js library
62 */
63Proj4js = {
64
65    /**
66     * Property: defaultDatum
67     * The datum to use when no others a specified
68     */
69    defaultDatum: 'WGS84',                  //default datum
70
71    /**
72    * Method: transform(source, dest, point)
73    * Transform a point coordinate from one map projection to another.  This is
74    * really the only public method you should need to use.
75    *
76    * Parameters:
77    * source - {Proj4js.Proj} source map projection for the transformation
78    * dest - {Proj4js.Proj} destination map projection for the transformation
79    * point - {Object} point to transform, may be geodetic (long, lat) or
80    *     projected Cartesian (x,y), but should always have x,y properties.
81    */
82    transform: function(source, dest, point) {
83        if (!source.readyToUse) {
84            this.reportError("Proj4js initialization for:"+source.srsCode+" not yet complete");
85            return point;
86        }
87        if (!dest.readyToUse) {
88            this.reportError("Proj4js initialization for:"+dest.srsCode+" not yet complete");
89            return point;
90        }
91       
92        // Workaround for Spherical Mercator
93        if ((source.srsProjNumber =="900913" && dest.datumCode != "WGS84" && !dest.datum_params) ||
94            (dest.srsProjNumber == "900913" && source.datumCode != "WGS84" && !source.datum_params)) {
95            var wgs84 = Proj4js.WGS84;
96            this.transform(source, wgs84, point);
97            source = wgs84;
98        }
99
100        // DGR, 2010/11/12
101        if (source.axis!="enu") {
102            this.adjust_axis(source,false,point);
103        }
104
105        // Transform source points to long/lat, if they aren't already.
106        if ( source.projName=="longlat") {
107            point.x *= Proj4js.common.D2R;  // convert degrees to radians
108            point.y *= Proj4js.common.D2R;
109        } else {
110            if (source.to_meter) {
111                point.x *= source.to_meter;
112                point.y *= source.to_meter;
113            }
114            source.inverse(point); // Convert Cartesian to longlat
115        }
116
117        // Adjust for the prime meridian if necessary
118        if (source.from_greenwich) { 
119            point.x += source.from_greenwich; 
120        }
121
122        // Convert datums if needed, and if possible.
123        point = this.datum_transform( source.datum, dest.datum, point );
124
125        // Adjust for the prime meridian if necessary
126        if (dest.from_greenwich) {
127            point.x -= dest.from_greenwich;
128        }
129
130        if( dest.projName=="longlat" ) {             
131            // convert radians to decimal degrees
132            point.x *= Proj4js.common.R2D;
133            point.y *= Proj4js.common.R2D;
134        } else  {               // else project
135            dest.forward(point);
136            if (dest.to_meter) {
137                point.x /= dest.to_meter;
138                point.y /= dest.to_meter;
139            }
140        }
141
142        // DGR, 2010/11/12
143        if (dest.axis!="enu") {
144            this.adjust_axis(dest,true,point);
145        }
146
147        return point;
148    }, // transform()
149
150    /** datum_transform()
151      source coordinate system definition,
152      destination coordinate system definition,
153      point to transform in geodetic coordinates (long, lat, height)
154    */
155    datum_transform : function( source, dest, point ) {
156
157      // Short cut if the datums are identical.
158      if( source.compare_datums( dest ) ) {
159          return point; // in this case, zero is sucess,
160                    // whereas cs_compare_datums returns 1 to indicate TRUE
161                    // confusing, should fix this
162      }
163
164      // Explicitly skip datum transform by setting 'datum=none' as parameter for either source or dest
165      if( source.datum_type == Proj4js.common.PJD_NODATUM
166          || dest.datum_type == Proj4js.common.PJD_NODATUM) {
167          return point;
168      }
169
170      // If this datum requires grid shifts, then apply it to geodetic coordinates.
171      if( source.datum_type == Proj4js.common.PJD_GRIDSHIFT )
172      {
173        alert("ERROR: Grid shift transformations are not implemented yet.");
174        /*
175          pj_apply_gridshift( pj_param(source.params,"snadgrids").s, 0,
176                              point_count, point_offset, x, y, z );
177          CHECK_RETURN;
178
179          src_a = SRS_WGS84_SEMIMAJOR;
180          src_es = 0.006694379990;
181        */
182      }
183
184      if( dest.datum_type == Proj4js.common.PJD_GRIDSHIFT )
185      {
186        alert("ERROR: Grid shift transformations are not implemented yet.");
187        /*
188          dst_a = ;
189          dst_es = 0.006694379990;
190        */
191      }
192
193      // Do we need to go through geocentric coordinates?
194      if( source.es != dest.es || source.a != dest.a
195          || source.datum_type == Proj4js.common.PJD_3PARAM
196          || source.datum_type == Proj4js.common.PJD_7PARAM
197          || dest.datum_type == Proj4js.common.PJD_3PARAM
198          || dest.datum_type == Proj4js.common.PJD_7PARAM)
199      {
200
201        // Convert to geocentric coordinates.
202        source.geodetic_to_geocentric( point );
203        // CHECK_RETURN;
204
205        // Convert between datums
206        if( source.datum_type == Proj4js.common.PJD_3PARAM || source.datum_type == Proj4js.common.PJD_7PARAM ) {
207          source.geocentric_to_wgs84(point);
208          // CHECK_RETURN;
209        }
210
211        if( dest.datum_type == Proj4js.common.PJD_3PARAM || dest.datum_type == Proj4js.common.PJD_7PARAM ) {
212          dest.geocentric_from_wgs84(point);
213          // CHECK_RETURN;
214        }
215
216        // Convert back to geodetic coordinates
217        dest.geocentric_to_geodetic( point );
218          // CHECK_RETURN;
219      }
220
221      // Apply grid shift to destination if required
222      if( dest.datum_type == Proj4js.common.PJD_GRIDSHIFT )
223      {
224        alert("ERROR: Grid shift transformations are not implemented yet.");
225        // pj_apply_gridshift( pj_param(dest.params,"snadgrids").s, 1, point);
226        // CHECK_RETURN;
227      }
228      return point;
229    }, // cs_datum_transform
230
231    /**
232     * Function: adjust_axis
233     * Normalize or de-normalized the x/y/z axes.  The normal form is "enu"
234     * (easting, northing, up).
235     * Parameters:
236     * crs {Proj4js.Proj} the coordinate reference system
237     * denorm {Boolean} when false, normalize
238     * point {Object} the coordinates to adjust
239     */
240    adjust_axis: function(crs, denorm, point) {
241        var xin= point.x, yin= point.y, zin= point.z || 0.0;
242        var v, t;
243        for (var i= 0; i<3; i++) {
244            if (denorm && i==2 && point.z===undefined) { continue; }
245                 if (i==0) { v= xin; t= 'x'; }
246            else if (i==1) { v= yin; t= 'y'; }
247            else           { v= zin; t= 'z'; }
248            switch(crs.axis[i]) {
249            case 'e':
250                point[t]= v;
251                break;
252            case 'w':
253                point[t]= -v;
254                break;
255            case 'n':
256                point[t]= v;
257                break;
258            case 's':
259                point[t]= -v;
260                break;
261            case 'u':
262                if (point[t]!==undefined) { point.z= v; }
263                break;
264            case 'd':
265                if (point[t]!==undefined) { point.z= -v; }
266                break;
267            default :
268                alert("ERROR: unknow axis ("+crs.axis[i]+") - check definition of "+src.projName);
269                return null;
270            }
271        }
272        return point;
273    },
274
275    /**
276     * Function: reportError
277     * An internal method to report errors back to user.
278     * Override this in applications to report error messages or throw exceptions.
279     */
280    reportError: function(msg) {
281      //console.log(msg);
282    },
283
284/**
285 *
286 * Title: Private Methods
287 * The following properties and methods are intended for internal use only.
288 *
289 * This is a minimal implementation of JavaScript inheritance methods so that
290 * Proj4js can be used as a stand-alone library.
291 * These are copies of the equivalent OpenLayers methods at v2.7
292 */
293 
294/**
295 * Function: extend
296 * Copy all properties of a source object to a destination object.  Modifies
297 *     the passed in destination object.  Any properties on the source object
298 *     that are set to undefined will not be (re)set on the destination object.
299 *
300 * Parameters:
301 * destination - {Object} The object that will be modified
302 * source - {Object} The object with properties to be set on the destination
303 *
304 * Returns:
305 * {Object} The destination object.
306 */
307    extend: function(destination, source) {
308      destination = destination || {};
309      if(source) {
310          for(var property in source) {
311              var value = source[property];
312              if(value !== undefined) {
313                  destination[property] = value;
314              }
315          }
316      }
317      return destination;
318    },
319
320/**
321 * Constructor: Class
322 * Base class used to construct all other classes. Includes support for
323 *     multiple inheritance.
324 * 
325 */
326    Class: function() {
327      var Class = function() {
328          this.initialize.apply(this, arguments);
329      };
330 
331      var extended = {};
332      var parent;
333      for(var i=0; i<arguments.length; ++i) {
334          if(typeof arguments[i] == "function") {
335              // get the prototype of the superclass
336              parent = arguments[i].prototype;
337          } else {
338              // in this case we're extending with the prototype
339              parent = arguments[i];
340          }
341          Proj4js.extend(extended, parent);
342      }
343      Class.prototype = extended;
344     
345      return Class;
346    },
347
348    /**
349     * Function: bind
350     * Bind a function to an object.  Method to easily create closures with
351     *     'this' altered.
352     *
353     * Parameters:
354     * func - {Function} Input function.
355     * object - {Object} The object to bind to the input function (as this).
356     *
357     * Returns:
358     * {Function} A closure with 'this' set to the passed in object.
359     */
360    bind: function(func, object) {
361        // create a reference to all arguments past the second one
362        var args = Array.prototype.slice.apply(arguments, [2]);
363        return function() {
364            // Push on any additional arguments from the actual function call.
365            // These will come after those sent to the bind call.
366            var newArgs = args.concat(
367                Array.prototype.slice.apply(arguments, [0])
368            );
369            return func.apply(object, newArgs);
370        };
371    },
372   
373/**
374 * The following properties and methods handle dynamic loading of JSON objects.
375 */
376 
377    /**
378     * Property: scriptName
379     * {String} The filename of this script without any path.
380     */
381    scriptName: "proj4js-combined.js",
382
383    /**
384     * Property: defsLookupService
385     * AJAX service to retreive projection definition parameters from
386     */
387    defsLookupService: 'http://spatialreference.org/ref',
388
389    /**
390     * Property: libPath
391     * internal: http server path to library code.
392     */
393    libPath: null,
394
395    /**
396     * Function: getScriptLocation
397     * Return the path to this script.
398     *
399     * Returns:
400     * Path to this script
401     */
402    getScriptLocation: function () {
403        if (this.libPath) return this.libPath;
404        var scriptName = this.scriptName;
405        var scriptNameLen = scriptName.length;
406
407        var scripts = document.getElementsByTagName('script');
408        for (var i = 0; i < scripts.length; i++) {
409            var src = scripts[i].getAttribute('src');
410            if (src) {
411                var index = src.lastIndexOf(scriptName);
412                // is it found, at the end of the URL?
413                if ((index > -1) && (index + scriptNameLen == src.length)) {
414                    this.libPath = src.slice(0, -scriptNameLen);
415                    break;
416                }
417            }
418        }
419        return this.libPath||"";
420    },
421
422    /**
423     * Function: loadScript
424     * Load a JS file from a URL into a <script> tag in the page.
425     *
426     * Parameters:
427     * url - {String} The URL containing the script to load
428     * onload - {Function} A method to be executed when the script loads successfully
429     * onfail - {Function} A method to be executed when there is an error loading the script
430     * loadCheck - {Function} A boolean method that checks to see if the script
431     *            has loaded.  Typically this just checks for the existance of
432     *            an object in the file just loaded.
433     */
434    loadScript: function(url, onload, onfail, loadCheck) {
435      var script = document.createElement('script');
436      script.defer = false;
437      script.type = "text/javascript";
438      script.id = url;
439      script.src = url;
440      script.onload = onload;
441      script.onerror = onfail;
442      script.loadCheck = loadCheck;
443      if (/MSIE/.test(navigator.userAgent)) {
444        script.onreadystatechange = this.checkReadyState;
445      }
446      document.getElementsByTagName('head')[0].appendChild(script);
447    },
448   
449    /**
450     * Function: checkReadyState
451     * IE workaround since there is no onerror handler.  Calls the user defined
452     * loadCheck method to determine if the script is loaded.
453     *
454     */
455    checkReadyState: function() {
456      if (this.readyState == 'loaded') {
457        if (!this.loadCheck()) {
458          this.onerror();
459        } else {
460          this.onload();
461        }
462      }
463    }
464};
465
466/**
467 * Class: Proj4js.Proj
468 *
469 * Proj objects provide transformation methods for point coordinates
470 * between geodetic latitude/longitude and a projected coordinate system.
471 * once they have been initialized with a projection code.
472 *
473 * Initialization of Proj objects is with a projection code, usually EPSG codes,
474 * which is the key that will be used with the Proj4js.defs array.
475 *
476 * The code passed in will be stripped of colons and converted to uppercase
477 * to locate projection definition files.
478 *
479 * A projection object has properties for units and title strings.
480 */
481Proj4js.Proj = Proj4js.Class({
482
483  /**
484   * Property: readyToUse
485   * Flag to indicate if initialization is complete for this Proj object
486   */
487  readyToUse: false,   
488 
489  /**
490   * Property: title
491   * The title to describe the projection
492   */
493  title: null, 
494 
495  /**
496   * Property: projName
497   * The projection class for this projection, e.g. lcc (lambert conformal conic,
498   * or merc for mercator).  These are exactly equivalent to their Proj4
499   * counterparts.
500   */
501  projName: null,
502  /**
503   * Property: units
504   * The units of the projection.  Values include 'm' and 'degrees'
505   */
506  units: null,
507  /**
508   * Property: datum
509   * The datum specified for the projection
510   */
511  datum: null,
512  /**
513   * Property: x0
514   * The x coordinate origin
515   */
516  x0: 0,
517  /**
518   * Property: y0
519   * The y coordinate origin
520   */
521  y0: 0,
522  /**
523   * Property: localCS
524   * Flag to indicate if the projection is a local one in which no transforms
525   * are required.
526   */
527  localCS: false,
528
529  /**
530  * Property: queue
531  * Buffer (FIFO) to hold callbacks waiting to be called when projection loaded.
532  */
533  queue: null,
534
535  /**
536  * Constructor: initialize
537  * Constructor for Proj4js.Proj objects
538  *
539  * Parameters:
540  * srsCode - a code for map projection definition parameters.  These are usually
541  * (but not always) EPSG codes.
542  */
543  initialize: function(srsCode, callback) {
544      this.srsCodeInput = srsCode;
545     
546      //Register callbacks prior to attempting to process definition
547      this.queue = [];
548      if( callback ){
549           this.queue.push( callback );
550      }
551     
552      //check to see if this is a WKT string
553      if ((srsCode.indexOf('GEOGCS') >= 0) ||
554          (srsCode.indexOf('GEOCCS') >= 0) ||
555          (srsCode.indexOf('PROJCS') >= 0) ||
556          (srsCode.indexOf('LOCAL_CS') >= 0)) {
557            this.parseWKT(srsCode);
558            this.deriveConstants();
559            this.loadProjCode(this.projName);
560            return;
561      }
562     
563      // DGR 2008-08-03 : support urn and url
564      if (srsCode.indexOf('urn:') == 0) {
565          //urn:ORIGINATOR:def:crs:CODESPACE:VERSION:ID
566          var urn = srsCode.split(':');
567          if ((urn[1] == 'ogc' || urn[1] =='x-ogc') &&
568              (urn[2] =='def') &&
569              (urn[3] =='crs')) {
570              srsCode = urn[4]+':'+urn[urn.length-1];
571          }
572      } else if (srsCode.indexOf('http://') == 0) {
573          //url#ID
574          var url = srsCode.split('#');
575          if (url[0].match(/epsg.org/)) {
576            // http://www.epsg.org/#
577            srsCode = 'EPSG:'+url[1];
578          } else if (url[0].match(/RIG.xml/)) {
579            //http://librairies.ign.fr/geoportail/resources/RIG.xml#
580            //http://interop.ign.fr/registers/ign/RIG.xml#
581            srsCode = 'IGNF:'+url[1];
582          }
583      }
584      this.srsCode = srsCode.toUpperCase();
585      if (this.srsCode.indexOf("EPSG") == 0) {
586          this.srsCode = this.srsCode;
587          this.srsAuth = 'epsg';
588          this.srsProjNumber = this.srsCode.substring(5);
589      // DGR 2007-11-20 : authority IGNF
590      } else if (this.srsCode.indexOf("IGNF") == 0) {
591          this.srsCode = this.srsCode;
592          this.srsAuth = 'IGNF';
593          this.srsProjNumber = this.srsCode.substring(5);
594      // DGR 2008-06-19 : pseudo-authority CRS for WMS
595      } else if (this.srsCode.indexOf("CRS") == 0) {
596          this.srsCode = this.srsCode;
597          this.srsAuth = 'CRS';
598          this.srsProjNumber = this.srsCode.substring(4);
599      } else {
600          this.srsAuth = '';
601          this.srsProjNumber = this.srsCode;
602      }
603     
604      this.loadProjDefinition();
605  },
606 
607/**
608 * Function: loadProjDefinition
609 *    Loads the coordinate system initialization string if required.
610 *    Note that dynamic loading happens asynchronously so an application must
611 *    wait for the readyToUse property is set to true.
612 *    To prevent dynamic loading, include the defs through a script tag in
613 *    your application.
614 *
615 */
616    loadProjDefinition: function() {
617      //check in memory
618      if (Proj4js.defs[this.srsCode]) {
619        this.defsLoaded();
620        return;
621      }
622
623      //else check for def on the server
624      var url = Proj4js.getScriptLocation() + 'defs/' + this.srsAuth.toUpperCase() + this.srsProjNumber + '.js';
625      Proj4js.loadScript(url, 
626                Proj4js.bind(this.defsLoaded, this),
627                Proj4js.bind(this.loadFromService, this),
628                Proj4js.bind(this.checkDefsLoaded, this) );
629    },
630
631/**
632 * Function: loadFromService
633 *    Creates the REST URL for loading the definition from a web service and
634 *    loads it.
635 *
636 */
637    loadFromService: function() {
638      //else load from web service
639      var url = Proj4js.defsLookupService +'/' + this.srsAuth +'/'+ this.srsProjNumber + '/proj4js/';
640      Proj4js.loadScript(url, 
641            Proj4js.bind(this.defsLoaded, this),
642            Proj4js.bind(this.defsFailed, this),
643            Proj4js.bind(this.checkDefsLoaded, this) );
644    },
645
646/**
647 * Function: defsLoaded
648 * Continues the Proj object initilization once the def file is loaded
649 *
650 */
651    defsLoaded: function() {
652      this.parseDefs();
653      this.loadProjCode(this.projName);
654    },
655   
656/**
657 * Function: checkDefsLoaded
658 *    This is the loadCheck method to see if the def object exists
659 *
660 */
661    checkDefsLoaded: function() {
662      if (Proj4js.defs[this.srsCode]) {
663        return true;
664      } else {
665        return false;
666      }
667    },
668
669 /**
670 * Function: defsFailed
671 *    Report an error in loading the defs file, but continue on using WGS84
672 *
673 */
674   defsFailed: function() {
675      Proj4js.reportError('failed to load projection definition for: '+this.srsCode);
676      Proj4js.defs[this.srsCode] = Proj4js.defs['WGS84'];  //set it to something so it can at least continue
677      this.defsLoaded();
678    },
679
680/**
681 * Function: loadProjCode
682 *    Loads projection class code dynamically if required.
683 *     Projection code may be included either through a script tag or in
684 *     a built version of proj4js
685 *
686 */
687    loadProjCode: function(projName) {
688      if (Proj4js.Proj[projName]) {
689        this.initTransforms();
690        return;
691      }
692
693      //the URL for the projection code
694      var url = Proj4js.getScriptLocation() + 'projCode/' + projName + '.js';
695      Proj4js.loadScript(url, 
696              Proj4js.bind(this.loadProjCodeSuccess, this, projName),
697              Proj4js.bind(this.loadProjCodeFailure, this, projName), 
698              Proj4js.bind(this.checkCodeLoaded, this, projName) );
699    },
700
701 /**
702 * Function: loadProjCodeSuccess
703 *    Loads any proj dependencies or continue on to final initialization.
704 *
705 */
706    loadProjCodeSuccess: function(projName) {
707      if (Proj4js.Proj[projName].dependsOn){
708        this.loadProjCode(Proj4js.Proj[projName].dependsOn);
709      } else {
710        this.initTransforms();
711      }
712    },
713
714 /**
715 * Function: defsFailed
716 *    Report an error in loading the proj file.  Initialization of the Proj
717 *    object has failed and the readyToUse flag will never be set.
718 *
719 */
720    loadProjCodeFailure: function(projName) {
721      Proj4js.reportError("failed to find projection file for: " + projName);
722      //TBD initialize with identity transforms so proj will still work?
723    },
724   
725/**
726 * Function: checkCodeLoaded
727 *    This is the loadCheck method to see if the projection code is loaded
728 *
729 */
730    checkCodeLoaded: function(projName) {
731      if (Proj4js.Proj[projName]) {
732        return true;
733      } else {
734        return false;
735      }
736    },
737
738/**
739 * Function: initTransforms
740 *    Finalize the initialization of the Proj object
741 *
742 */
743    initTransforms: function() {
744      Proj4js.extend(this, Proj4js.Proj[this.projName]);
745      this.init();
746      this.readyToUse = true;
747      if( this.queue ) {
748        var item;
749        while( (item = this.queue.shift()) ) {
750          item.call( this, this );
751        }
752      }
753  },
754
755/**
756 * Function: parseWKT
757 * Parses a WKT string to get initialization parameters
758 *
759 */
760 wktRE: /^(\w+)\[(.*)\]$/,
761 parseWKT: function(wkt) {
762    var wktMatch = wkt.match(this.wktRE);
763    if (!wktMatch) return;
764    var wktObject = wktMatch[1];
765    var wktContent = wktMatch[2];
766    var wktTemp = wktContent.split(",");
767    var wktName;
768    if (wktObject.toUpperCase() == "TOWGS84") {
769      wktName = wktObject;  //no name supplied for the TOWGS84 array
770    } else {
771      wktName = wktTemp.shift();
772    }
773    wktName = wktName.replace(/^\"/,"");
774    wktName = wktName.replace(/\"$/,"");
775   
776    /*
777    wktContent = wktTemp.join(",");
778    var wktArray = wktContent.split("],");
779    for (var i=0; i<wktArray.length-1; ++i) {
780      wktArray[i] += "]";
781    }
782    */
783   
784    var wktArray = new Array();
785    var bkCount = 0;
786    var obj = "";
787    for (var i=0; i<wktTemp.length; ++i) {
788      var token = wktTemp[i];
789      for (var j=0; j<token.length; ++j) {
790        if (token.charAt(j) == "[") ++bkCount;
791        if (token.charAt(j) == "]") --bkCount;
792      }
793      obj += token;
794      if (bkCount === 0) {
795        wktArray.push(obj);
796        obj = "";
797      } else {
798        obj += ",";
799      }
800    }
801   
802    //do something based on the type of the wktObject being parsed
803    //add in variations in the spelling as required
804    switch (wktObject) {
805      case 'LOCAL_CS':
806        this.projName = 'identity'
807        this.localCS = true;
808        this.srsCode = wktName;
809        break;
810      case 'GEOGCS':
811        this.projName = 'longlat'
812        this.geocsCode = wktName;
813        if (!this.srsCode) this.srsCode = wktName;
814        break;
815      case 'PROJCS':
816        this.srsCode = wktName;
817        break;
818      case 'GEOCCS':
819        break;
820      case 'PROJECTION':
821        this.projName = Proj4js.wktProjections[wktName]
822        break;
823      case 'DATUM':
824        this.datumName = wktName;
825        break;
826      case 'LOCAL_DATUM':
827        this.datumCode = 'none';
828        break;
829      case 'SPHEROID':
830        this.ellps = wktName;
831        this.a = parseFloat(wktArray.shift());
832        this.rf = parseFloat(wktArray.shift());
833        break;
834      case 'PRIMEM':
835        this.from_greenwich = parseFloat(wktArray.shift()); //to radians?
836        break;
837      case 'UNIT':
838        this.units = wktName;
839        this.unitsPerMeter = parseFloat(wktArray.shift());
840        break;
841      case 'PARAMETER':
842        var name = wktName.toLowerCase();
843        var value = parseFloat(wktArray.shift());
844        //there may be many variations on the wktName values, add in case
845        //statements as required
846        switch (name) {
847          case 'false_easting':
848            this.x0 = value;
849            break;
850          case 'false_northing':
851            this.y0 = value;
852            break;
853          case 'scale_factor':
854            this.k0 = value;
855            break;
856          case 'central_meridian':
857            this.long0 = value*Proj4js.common.D2R;
858            break;
859          case 'latitude_of_origin':
860            this.lat0 = value*Proj4js.common.D2R;
861            break;
862          case 'more_here':
863            break;
864          default:
865            break;
866        }
867        break;
868      case 'TOWGS84':
869        this.datum_params = wktArray;
870        break;
871      //DGR 2010-11-12: AXIS
872      case 'AXIS':
873        var name= wktName.toLowerCase();
874        var value= wktArray.shift();
875        switch (value) {
876          case 'EAST' : value= 'e'; break;
877          case 'WEST' : value= 'w'; break;
878          case 'NORTH': value= 'n'; break;
879          case 'SOUTH': value= 's'; break;
880          case 'UP'   : value= 'u'; break;
881          case 'DOWN' : value= 'd'; break;
882          case 'OTHER':
883          default     : value= ' '; break;//FIXME
884        }
885        if (!this.axis) { this.axis= "enu"; }
886        switch(name) {
887          case 'X': this.axis=                         value + this.axis.substr(1,2); break;
888          case 'Y': this.axis= this.axis.substr(0,1) + value + this.axis.substr(2,1); break;
889          case 'Z': this.axis= this.axis.substr(0,2) + value                        ; break;
890          default : break;
891        }
892      case 'MORE_HERE':
893        break;
894      default:
895        break;
896    }
897    for (var i=0; i<wktArray.length; ++i) {
898      this.parseWKT(wktArray[i]);
899    }
900 },
901
902/**
903 * Function: parseDefs
904 * Parses the PROJ.4 initialization string and sets the associated properties.
905 *
906 */
907  parseDefs: function() {
908      this.defData = Proj4js.defs[this.srsCode];
909      var paramName, paramVal;
910      if (!this.defData) {
911        return;
912      }
913      var paramArray=this.defData.split("+");
914
915      for (var prop=0; prop<paramArray.length; prop++) {
916          var property = paramArray[prop].split("=");
917          paramName = property[0].toLowerCase();
918          paramVal = property[1];
919
920          switch (paramName.replace(/\s/gi,"")) {  // trim out spaces
921              case "": break;   // throw away nameless parameter
922              case "title":  this.title = paramVal; break;
923              case "proj":   this.projName =  paramVal.replace(/\s/gi,""); break;
924              case "units":  this.units = paramVal.replace(/\s/gi,""); break;
925              case "datum":  this.datumCode = paramVal.replace(/\s/gi,""); break;
926              case "nadgrids": this.nagrids = paramVal.replace(/\s/gi,""); break;
927              case "ellps":  this.ellps = paramVal.replace(/\s/gi,""); break;
928              case "a":      this.a =  parseFloat(paramVal); break;  // semi-major radius
929              case "b":      this.b =  parseFloat(paramVal); break;  // semi-minor radius
930              // DGR 2007-11-20
931              case "rf":     this.rf = parseFloat(paramVal); break; // inverse flattening rf= a/(a-b)
932              case "lat_0":  this.lat0 = paramVal*Proj4js.common.D2R; break;        // phi0, central latitude
933              case "lat_1":  this.lat1 = paramVal*Proj4js.common.D2R; break;        //standard parallel 1
934              case "lat_2":  this.lat2 = paramVal*Proj4js.common.D2R; break;        //standard parallel 2
935              case "lat_ts": this.lat_ts = paramVal*Proj4js.common.D2R; break;      // used in merc and eqc
936              case "lon_0":  this.long0 = paramVal*Proj4js.common.D2R; break;       // lam0, central longitude
937              case "alpha":  this.alpha =  parseFloat(paramVal)*Proj4js.common.D2R; break;  //for somerc projection
938              case "lonc":   this.longc = paramVal*Proj4js.common.D2R; break;       //for somerc projection
939              case "x_0":    this.x0 = parseFloat(paramVal); break;  // false easting
940              case "y_0":    this.y0 = parseFloat(paramVal); break;  // false northing
941              case "k_0":    this.k0 = parseFloat(paramVal); break;  // projection scale factor
942              case "k":      this.k0 = parseFloat(paramVal); break;  // both forms returned
943              case "r_a":    this.R_A = true; break;                 // sphere--area of ellipsoid
944              case "zone":   this.zone = parseInt(paramVal); break;  // UTM Zone
945              case "south":   this.utmSouth = true; break;  // UTM north/south
946              case "towgs84":this.datum_params = paramVal.split(","); break;
947              case "to_meter": this.to_meter = parseFloat(paramVal); break; // cartesian scaling
948              case "from_greenwich": this.from_greenwich = paramVal*Proj4js.common.D2R; break;
949              // DGR 2008-07-09 : if pm is not a well-known prime meridian take
950              // the value instead of 0.0, then convert to radians
951              case "pm":     paramVal = paramVal.replace(/\s/gi,"");
952                             this.from_greenwich = Proj4js.PrimeMeridian[paramVal] ?
953                                Proj4js.PrimeMeridian[paramVal] : parseFloat(paramVal);
954                             this.from_greenwich *= Proj4js.common.D2R; 
955                             break;
956              // DGR 2010-11-12: axis
957              case "axis":   paramVal = paramVal.replace(/\s/gi,"");
958                             var legalAxis= "ewnsud";
959                             if (paramVal.length==3 &&
960                                 legalAxis.indexOf(paramVal.substr(0,1))!=-1 &&
961                                 legalAxis.indexOf(paramVal.substr(1,1))!=-1 &&
962                                 legalAxis.indexOf(paramVal.substr(2,1))!=-1) {
963                                this.axis= paramVal;
964                             } //FIXME: be silent ?
965                             break
966              case "no_defs": break; 
967              default: //alert("Unrecognized parameter: " + paramName);
968          } // switch()
969      } // for paramArray
970      this.deriveConstants();
971  },
972
973/**
974 * Function: deriveConstants
975 * Sets several derived constant values and initialization of datum and ellipse
976 *     parameters.
977 *
978 */
979  deriveConstants: function() {
980      if (this.nagrids == '@null') this.datumCode = 'none';
981      if (this.datumCode && this.datumCode != 'none') {
982        var datumDef = Proj4js.Datum[this.datumCode];
983        if (datumDef) {
984          this.datum_params = datumDef.towgs84 ? datumDef.towgs84.split(',') : null;
985          this.ellps = datumDef.ellipse;
986          this.datumName = datumDef.datumName ? datumDef.datumName : this.datumCode;
987        }
988      }
989      if (!this.a) {    // do we have an ellipsoid?
990          var ellipse = Proj4js.Ellipsoid[this.ellps] ? Proj4js.Ellipsoid[this.ellps] : Proj4js.Ellipsoid['WGS84'];
991          Proj4js.extend(this, ellipse);
992      }
993      if (this.rf && !this.b) this.b = (1.0 - 1.0/this.rf) * this.a;
994      if (Math.abs(this.a - this.b)<Proj4js.common.EPSLN) {
995        this.sphere = true;
996        this.b= this.a;
997      }
998      this.a2 = this.a * this.a;          // used in geocentric
999      this.b2 = this.b * this.b;          // used in geocentric
1000      this.es = (this.a2-this.b2)/this.a2;  // e ^ 2
1001      this.e = Math.sqrt(this.es);        // eccentricity
1002      if (this.R_A) {
1003        this.a *= 1. - this.es * (Proj4js.common.SIXTH + this.es * (Proj4js.common.RA4 + this.es * Proj4js.common.RA6));
1004        this.a2 = this.a * this.a;
1005        this.b2 = this.b * this.b;
1006        this.es = 0.;
1007      }
1008      this.ep2=(this.a2-this.b2)/this.b2; // used in geocentric
1009      if (!this.k0) this.k0 = 1.0;    //default value
1010      //DGR 2010-11-12: axis
1011      if (!this.axis) { this.axis= "enu"; }
1012
1013      this.datum = new Proj4js.datum(this);
1014  }
1015});
1016
1017Proj4js.Proj.longlat = {
1018  init: function() {
1019    //no-op for longlat
1020  },
1021  forward: function(pt) {
1022    //identity transform
1023    return pt;
1024  },
1025  inverse: function(pt) {
1026    //identity transform
1027    return pt;
1028  }
1029};
1030Proj4js.Proj.identity = Proj4js.Proj.longlat;
1031
1032/**
1033  Proj4js.defs is a collection of coordinate system definition objects in the
1034  PROJ.4 command line format.
1035  Generally a def is added by means of a separate .js file for example:
1036
1037    <SCRIPT type="text/javascript" src="defs/EPSG26912.js"></SCRIPT>
1038
1039  def is a CS definition in PROJ.4 WKT format, for example:
1040    +proj="tmerc"   //longlat, etc.
1041    +a=majorRadius
1042    +b=minorRadius
1043    +lat0=somenumber
1044    +long=somenumber
1045*/
1046Proj4js.defs = {
1047  // These are so widely used, we'll go ahead and throw them in
1048  // without requiring a separate .js file
1049  'WGS84': "+title=long/lat:WGS84 +proj=longlat +ellps=WGS84 +datum=WGS84 +units=degrees",
1050  'EPSG:4326': "+title=long/lat:WGS84 +proj=longlat +a=6378137.0 +b=6356752.31424518 +ellps=WGS84 +datum=WGS84 +units=degrees",
1051  'EPSG:4269': "+title=long/lat:NAD83 +proj=longlat +a=6378137.0 +b=6356752.31414036 +ellps=GRS80 +datum=NAD83 +units=degrees",
1052  'EPSG:3785': "+title= Google Mercator +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"
1053};
1054Proj4js.defs['GOOGLE'] = Proj4js.defs['EPSG:3785'];
1055Proj4js.defs['EPSG:900913'] = Proj4js.defs['EPSG:3785'];
1056Proj4js.defs['EPSG:102113'] = Proj4js.defs['EPSG:3785'];
1057
1058Proj4js.common = {
1059  PI : 3.141592653589793238, //Math.PI,
1060  HALF_PI : 1.570796326794896619, //Math.PI*0.5,
1061  TWO_PI : 6.283185307179586477, //Math.PI*2,
1062  FORTPI : 0.78539816339744833,
1063  R2D : 57.29577951308232088,
1064  D2R : 0.01745329251994329577,
1065  SEC_TO_RAD : 4.84813681109535993589914102357e-6, /* SEC_TO_RAD = Pi/180/3600 */
1066  EPSLN : 1.0e-10,
1067  MAX_ITER : 20,
1068  // following constants from geocent.c
1069  COS_67P5 : 0.38268343236508977,  /* cosine of 67.5 degrees */
1070  AD_C : 1.0026000,                /* Toms region 1 constant */
1071
1072  /* datum_type values */
1073  PJD_UNKNOWN  : 0,
1074  PJD_3PARAM   : 1,
1075  PJD_7PARAM   : 2,
1076  PJD_GRIDSHIFT: 3,
1077  PJD_WGS84    : 4,   // WGS84 or equivalent
1078  PJD_NODATUM  : 5,   // WGS84 or equivalent
1079  SRS_WGS84_SEMIMAJOR : 6378137.0,  // only used in grid shift transforms
1080
1081  // ellipoid pj_set_ell.c
1082  SIXTH : .1666666666666666667, /* 1/6 */
1083  RA4   : .04722222222222222222, /* 17/360 */
1084  RA6   : .02215608465608465608, /* 67/3024 */
1085  RV4   : .06944444444444444444, /* 5/72 */
1086  RV6   : .04243827160493827160, /* 55/1296 */
1087
1088// Function to compute the constant small m which is the radius of
1089//   a parallel of latitude, phi, divided by the semimajor axis.
1090// -----------------------------------------------------------------
1091  msfnz : function(eccent, sinphi, cosphi) {
1092      var con = eccent * sinphi;
1093      return cosphi/(Math.sqrt(1.0 - con * con));
1094  },
1095
1096// Function to compute the constant small t for use in the forward
1097//   computations in the Lambert Conformal Conic and the Polar
1098//   Stereographic projections.
1099// -----------------------------------------------------------------
1100  tsfnz : function(eccent, phi, sinphi) {
1101    var con = eccent * sinphi;
1102    var com = .5 * eccent;
1103    con = Math.pow(((1.0 - con) / (1.0 + con)), com);
1104    return (Math.tan(.5 * (this.HALF_PI - phi))/con);
1105  },
1106
1107// Function to compute the latitude angle, phi2, for the inverse of the
1108//   Lambert Conformal Conic and Polar Stereographic projections.
1109// ----------------------------------------------------------------
1110  phi2z : function(eccent, ts) {
1111    var eccnth = .5 * eccent;
1112    var con, dphi;
1113    var phi = this.HALF_PI - 2 * Math.atan(ts);
1114    for (var i = 0; i <= 15; i++) {
1115      con = eccent * Math.sin(phi);
1116      dphi = this.HALF_PI - 2 * Math.atan(ts *(Math.pow(((1.0 - con)/(1.0 + con)),eccnth))) - phi;
1117      phi += dphi;
1118      if (Math.abs(dphi) <= .0000000001) return phi;
1119    }
1120    alert("phi2z has NoConvergence");
1121    return (-9999);
1122  },
1123
1124/* Function to compute constant small q which is the radius of a
1125   parallel of latitude, phi, divided by the semimajor axis.
1126------------------------------------------------------------*/
1127  qsfnz : function(eccent,sinphi) {
1128    var con;
1129    if (eccent > 1.0e-7) {
1130      con = eccent * sinphi;
1131      return (( 1.0- eccent * eccent) * (sinphi /(1.0 - con * con) - (.5/eccent)*Math.log((1.0 - con)/(1.0 + con))));
1132    } else {
1133      return(2.0 * sinphi);
1134    }
1135  },
1136
1137/* Function to eliminate roundoff errors in asin
1138----------------------------------------------*/
1139  asinz : function(x) {
1140    if (Math.abs(x)>1.0) {
1141      x=(x>1.0)?1.0:-1.0;
1142    }
1143    return Math.asin(x);
1144  },
1145
1146// following functions from gctpc cproj.c for transverse mercator projections
1147  e0fn : function(x) {return(1.0-0.25*x*(1.0+x/16.0*(3.0+1.25*x)));},
1148  e1fn : function(x) {return(0.375*x*(1.0+0.25*x*(1.0+0.46875*x)));},
1149  e2fn : function(x) {return(0.05859375*x*x*(1.0+0.75*x));},
1150  e3fn : function(x) {return(x*x*x*(35.0/3072.0));},
1151  mlfn : function(e0,e1,e2,e3,phi) {return(e0*phi-e1*Math.sin(2.0*phi)+e2*Math.sin(4.0*phi)-e3*Math.sin(6.0*phi));},
1152
1153  srat : function(esinp, exp) {
1154    return(Math.pow((1.0-esinp)/(1.0+esinp), exp));
1155  },
1156
1157// Function to return the sign of an argument
1158  sign : function(x) { if (x < 0.0) return(-1); else return(1);},
1159
1160// Function to adjust longitude to -180 to 180; input in radians
1161  adjust_lon : function(x) {
1162    x = (Math.abs(x) < this.PI) ? x: (x - (this.sign(x)*this.TWO_PI) );
1163    return x;
1164  },
1165
1166// IGNF - DGR : algorithms used by IGN France
1167
1168// Function to adjust latitude to -90 to 90; input in radians
1169  adjust_lat : function(x) {
1170    x= (Math.abs(x) < this.HALF_PI) ? x: (x - (this.sign(x)*this.PI) );
1171    return x;
1172  },
1173
1174// Latitude Isometrique - close to tsfnz ...
1175  latiso : function(eccent, phi, sinphi) {
1176    if (Math.abs(phi) > this.HALF_PI) return +Number.NaN;
1177    if (phi==this.HALF_PI) return Number.POSITIVE_INFINITY;
1178    if (phi==-1.0*this.HALF_PI) return -1.0*Number.POSITIVE_INFINITY;
1179
1180    var con= eccent*sinphi;
1181    return Math.log(Math.tan((this.HALF_PI+phi)/2.0))+eccent*Math.log((1.0-con)/(1.0+con))/2.0;
1182  },
1183
1184  fL : function(x,L) {
1185    return 2.0*Math.atan(x*Math.exp(L)) - this.HALF_PI;
1186  },
1187
1188// Inverse Latitude Isometrique - close to ph2z
1189  invlatiso : function(eccent, ts) {
1190    var phi= this.fL(1.0,ts);
1191    var Iphi= 0.0;
1192    var con= 0.0;
1193    do {
1194      Iphi= phi;
1195      con= eccent*Math.sin(Iphi);
1196      phi= this.fL(Math.exp(eccent*Math.log((1.0+con)/(1.0-con))/2.0),ts)
1197    } while (Math.abs(phi-Iphi)>1.0e-12);
1198    return phi;
1199  },
1200
1201// Needed for Gauss Schreiber
1202// Original:  Denis Makarov (info@binarythings.com)
1203// Web Site:  http://www.binarythings.com
1204  sinh : function(x)
1205  {
1206    var r= Math.exp(x);
1207    r= (r-1.0/r)/2.0;
1208    return r;
1209  },
1210
1211  cosh : function(x)
1212  {
1213    var r= Math.exp(x);
1214    r= (r+1.0/r)/2.0;
1215    return r;
1216  },
1217
1218  tanh : function(x)
1219  {
1220    var r= Math.exp(x);
1221    r= (r-1.0/r)/(r+1.0/r);
1222    return r;
1223  },
1224
1225  asinh : function(x)
1226  {
1227    var s= (x>= 0? 1.0:-1.0);
1228    return s*(Math.log( Math.abs(x) + Math.sqrt(x*x+1.0) ));
1229  },
1230
1231  acosh : function(x)
1232  {
1233    return 2.0*Math.log(Math.sqrt((x+1.0)/2.0) + Math.sqrt((x-1.0)/2.0));
1234  },
1235
1236  atanh : function(x)
1237  {
1238    return Math.log((x-1.0)/(x+1.0))/2.0;
1239  },
1240
1241// Grande Normale
1242  gN : function(a,e,sinphi)
1243  {
1244    var temp= e*sinphi;
1245    return a/Math.sqrt(1.0 - temp*temp);
1246  }
1247
1248};
1249
1250/** datum object
1251*/
1252Proj4js.datum = Proj4js.Class({
1253
1254  initialize : function(proj) {
1255    this.datum_type = Proj4js.common.PJD_WGS84;   //default setting
1256    if (proj.datumCode && proj.datumCode == 'none') {
1257      this.datum_type = Proj4js.common.PJD_NODATUM;
1258    }
1259    if (proj && proj.datum_params) {
1260      for (var i=0; i<proj.datum_params.length; i++) {
1261        proj.datum_params[i]=parseFloat(proj.datum_params[i]);
1262      }
1263      if (proj.datum_params[0] != 0 || proj.datum_params[1] != 0 || proj.datum_params[2] != 0 ) {
1264        this.datum_type = Proj4js.common.PJD_3PARAM;
1265      }
1266      if (proj.datum_params.length > 3) {
1267        if (proj.datum_params[3] != 0 || proj.datum_params[4] != 0 ||
1268            proj.datum_params[5] != 0 || proj.datum_params[6] != 0 ) {
1269          this.datum_type = Proj4js.common.PJD_7PARAM;
1270          proj.datum_params[3] *= Proj4js.common.SEC_TO_RAD;
1271          proj.datum_params[4] *= Proj4js.common.SEC_TO_RAD;
1272          proj.datum_params[5] *= Proj4js.common.SEC_TO_RAD;
1273          proj.datum_params[6] = (proj.datum_params[6]/1000000.0) + 1.0;
1274        }
1275      }
1276    }
1277    if (proj) {
1278      this.a = proj.a;    //datum object also uses these values
1279      this.b = proj.b;
1280      this.es = proj.es;
1281      this.ep2 = proj.ep2;
1282      this.datum_params = proj.datum_params;
1283    }
1284  },
1285
1286  /****************************************************************/
1287  // cs_compare_datums()
1288  //   Returns 1 (TRUE) if the two datums match, otherwise 0 (FALSE).
1289  compare_datums : function( dest ) {
1290    if( this.datum_type != dest.datum_type ) {
1291      return false; // false, datums are not equal
1292    } else if( this.a != dest.a || Math.abs(this.es-dest.es) > 0.000000000050 ) {
1293      // the tolerence for es is to ensure that GRS80 and WGS84
1294      // are considered identical
1295      return false;
1296    } else if( this.datum_type == Proj4js.common.PJD_3PARAM ) {
1297      return (this.datum_params[0] == dest.datum_params[0]
1298              && this.datum_params[1] == dest.datum_params[1]
1299              && this.datum_params[2] == dest.datum_params[2]);
1300    } else if( this.datum_type == Proj4js.common.PJD_7PARAM ) {
1301      return (this.datum_params[0] == dest.datum_params[0]
1302              && this.datum_params[1] == dest.datum_params[1]
1303              && this.datum_params[2] == dest.datum_params[2]
1304              && this.datum_params[3] == dest.datum_params[3]
1305              && this.datum_params[4] == dest.datum_params[4]
1306              && this.datum_params[5] == dest.datum_params[5]
1307              && this.datum_params[6] == dest.datum_params[6]);
1308    } else if( this.datum_type == Proj4js.common.PJD_GRIDSHIFT ) {
1309      return strcmp( pj_param(this.params,"snadgrids").s,
1310                     pj_param(dest.params,"snadgrids").s ) == 0;
1311    } else {
1312      return true; // datums are equal
1313    }
1314  }, // cs_compare_datums()
1315
1316  /*
1317   * The function Convert_Geodetic_To_Geocentric converts geodetic coordinates
1318   * (latitude, longitude, and height) to geocentric coordinates (X, Y, Z),
1319   * according to the current ellipsoid parameters.
1320   *
1321   *    Latitude  : Geodetic latitude in radians                     (input)
1322   *    Longitude : Geodetic longitude in radians                    (input)
1323   *    Height    : Geodetic height, in meters                       (input)
1324   *    X         : Calculated Geocentric X coordinate, in meters    (output)
1325   *    Y         : Calculated Geocentric Y coordinate, in meters    (output)
1326   *    Z         : Calculated Geocentric Z coordinate, in meters    (output)
1327   *
1328   */
1329  geodetic_to_geocentric : function(p) {
1330    var Longitude = p.x;
1331    var Latitude = p.y;
1332    var Height = p.z ? p.z : 0;   //Z value not always supplied
1333    var X;  // output
1334    var Y;
1335    var Z;
1336
1337    var Error_Code=0;  //  GEOCENT_NO_ERROR;
1338    var Rn;            /*  Earth radius at location  */
1339    var Sin_Lat;       /*  Math.sin(Latitude)  */
1340    var Sin2_Lat;      /*  Square of Math.sin(Latitude)  */
1341    var Cos_Lat;       /*  Math.cos(Latitude)  */
1342
1343    /*
1344    ** Don't blow up if Latitude is just a little out of the value
1345    ** range as it may just be a rounding issue.  Also removed longitude
1346    ** test, it should be wrapped by Math.cos() and Math.sin().  NFW for PROJ.4, Sep/2001.
1347    */
1348    if( Latitude < -Proj4js.common.HALF_PI && Latitude > -1.001 * Proj4js.common.HALF_PI ) {
1349        Latitude = -Proj4js.common.HALF_PI;
1350    } else if( Latitude > Proj4js.common.HALF_PI && Latitude < 1.001 * Proj4js.common.HALF_PI ) {
1351        Latitude = Proj4js.common.HALF_PI;
1352    } else if ((Latitude < -Proj4js.common.HALF_PI) || (Latitude > Proj4js.common.HALF_PI)) {
1353      /* Latitude out of range */
1354      Proj4js.reportError('geocent:lat out of range:'+Latitude);
1355      return null;
1356    }
1357
1358    if (Longitude > Proj4js.common.PI) Longitude -= (2*Proj4js.common.PI);
1359    Sin_Lat = Math.sin(Latitude);
1360    Cos_Lat = Math.cos(Latitude);
1361    Sin2_Lat = Sin_Lat * Sin_Lat;
1362    Rn = this.a / (Math.sqrt(1.0e0 - this.es * Sin2_Lat));
1363    X = (Rn + Height) * Cos_Lat * Math.cos(Longitude);
1364    Y = (Rn + Height) * Cos_Lat * Math.sin(Longitude);
1365    Z = ((Rn * (1 - this.es)) + Height) * Sin_Lat;
1366
1367    p.x = X;
1368    p.y = Y;
1369    p.z = Z;
1370    return Error_Code;
1371  }, // cs_geodetic_to_geocentric()
1372
1373
1374  geocentric_to_geodetic : function (p) {
1375/* local defintions and variables */
1376/* end-criterium of loop, accuracy of sin(Latitude) */
1377var genau = 1.E-12;
1378var genau2 = (genau*genau);
1379var maxiter = 30;
1380
1381    var P;        /* distance between semi-minor axis and location */
1382    var RR;       /* distance between center and location */
1383    var CT;       /* sin of geocentric latitude */
1384    var ST;       /* cos of geocentric latitude */
1385    var RX;
1386    var RK;
1387    var RN;       /* Earth radius at location */
1388    var CPHI0;    /* cos of start or old geodetic latitude in iterations */
1389    var SPHI0;    /* sin of start or old geodetic latitude in iterations */
1390    var CPHI;     /* cos of searched geodetic latitude */
1391    var SPHI;     /* sin of searched geodetic latitude */
1392    var SDPHI;    /* end-criterium: addition-theorem of sin(Latitude(iter)-Latitude(iter-1)) */
1393    var At_Pole;     /* indicates location is in polar region */
1394    var iter;        /* # of continous iteration, max. 30 is always enough (s.a.) */
1395
1396    var X = p.x;
1397    var Y = p.y;
1398    var Z = p.z ? p.z : 0.0;   //Z value not always supplied
1399    var Longitude;
1400    var Latitude;
1401    var Height;
1402
1403    At_Pole = false;
1404    P = Math.sqrt(X*X+Y*Y);
1405    RR = Math.sqrt(X*X+Y*Y+Z*Z);
1406
1407/*      special cases for latitude and longitude */
1408    if (P/this.a < genau) {
1409
1410/*  special case, if P=0. (X=0., Y=0.) */
1411        At_Pole = true;
1412        Longitude = 0.0;
1413
1414/*  if (X,Y,Z)=(0.,0.,0.) then Height becomes semi-minor axis
1415 *  of ellipsoid (=center of mass), Latitude becomes PI/2 */
1416        if (RR/this.a < genau) {
1417            Latitude = Proj4js.common.HALF_PI;
1418            Height   = -this.b;
1419            return;
1420        }
1421    } else {
1422/*  ellipsoidal (geodetic) longitude
1423 *  interval: -PI < Longitude <= +PI */
1424        Longitude=Math.atan2(Y,X);
1425    }
1426
1427/* --------------------------------------------------------------
1428 * Following iterative algorithm was developped by
1429 * "Institut fï¿œr Erdmessung", University of Hannover, July 1988.
1430 * Internet: www.ife.uni-hannover.de
1431 * Iterative computation of CPHI,SPHI and Height.
1432 * Iteration of CPHI and SPHI to 10**-12 radian resp.
1433 * 2*10**-7 arcsec.
1434 * --------------------------------------------------------------
1435 */
1436    CT = Z/RR;
1437    ST = P/RR;
1438    RX = 1.0/Math.sqrt(1.0-this.es*(2.0-this.es)*ST*ST);
1439    CPHI0 = ST*(1.0-this.es)*RX;
1440    SPHI0 = CT*RX;
1441    iter = 0;
1442
1443/* loop to find sin(Latitude) resp. Latitude
1444 * until |sin(Latitude(iter)-Latitude(iter-1))| < genau */
1445    do
1446    {
1447        iter++;
1448        RN = this.a/Math.sqrt(1.0-this.es*SPHI0*SPHI0);
1449
1450/*  ellipsoidal (geodetic) height */
1451        Height = P*CPHI0+Z*SPHI0-RN*(1.0-this.es*SPHI0*SPHI0);
1452
1453        RK = this.es*RN/(RN+Height);
1454        RX = 1.0/Math.sqrt(1.0-RK*(2.0-RK)*ST*ST);
1455        CPHI = ST*(1.0-RK)*RX;
1456        SPHI = CT*RX;
1457        SDPHI = SPHI*CPHI0-CPHI*SPHI0;
1458        CPHI0 = CPHI;
1459        SPHI0 = SPHI;
1460    }
1461    while (SDPHI*SDPHI > genau2 && iter < maxiter);
1462
1463/*      ellipsoidal (geodetic) latitude */
1464    Latitude=Math.atan(SPHI/Math.abs(CPHI));
1465
1466    p.x = Longitude;
1467    p.y = Latitude;
1468    p.z = Height;
1469    return p;
1470  }, // cs_geocentric_to_geodetic()
1471
1472  /** Convert_Geocentric_To_Geodetic
1473   * The method used here is derived from 'An Improved Algorithm for
1474   * Geocentric to Geodetic Coordinate Conversion', by Ralph Toms, Feb 1996
1475   */
1476  geocentric_to_geodetic_noniter : function (p) {
1477    var X = p.x;
1478    var Y = p.y;
1479    var Z = p.z ? p.z : 0;   //Z value not always supplied
1480    var Longitude;
1481    var Latitude;
1482    var Height;
1483
1484    var W;        /* distance from Z axis */
1485    var W2;       /* square of distance from Z axis */
1486    var T0;       /* initial estimate of vertical component */
1487    var T1;       /* corrected estimate of vertical component */
1488    var S0;       /* initial estimate of horizontal component */
1489    var S1;       /* corrected estimate of horizontal component */
1490    var Sin_B0;   /* Math.sin(B0), B0 is estimate of Bowring aux variable */
1491    var Sin3_B0;  /* cube of Math.sin(B0) */
1492    var Cos_B0;   /* Math.cos(B0) */
1493    var Sin_p1;   /* Math.sin(phi1), phi1 is estimated latitude */
1494    var Cos_p1;   /* Math.cos(phi1) */
1495    var Rn;       /* Earth radius at location */
1496    var Sum;      /* numerator of Math.cos(phi1) */
1497    var At_Pole;  /* indicates location is in polar region */
1498
1499    X = parseFloat(X);  // cast from string to float
1500    Y = parseFloat(Y);
1501    Z = parseFloat(Z);
1502
1503    At_Pole = false;
1504    if (X != 0.0)
1505    {
1506        Longitude = Math.atan2(Y,X);
1507    }
1508    else
1509    {
1510        if (Y > 0)
1511        {
1512            Longitude = Proj4js.common.HALF_PI;
1513        }
1514        else if (Y < 0)
1515        {
1516            Longitude = -Proj4js.common.HALF_PI;
1517        }
1518        else
1519        {
1520            At_Pole = true;
1521            Longitude = 0.0;
1522            if (Z > 0.0)
1523            {  /* north pole */
1524                Latitude = Proj4js.common.HALF_PI;
1525            }
1526            else if (Z < 0.0)
1527            {  /* south pole */
1528                Latitude = -Proj4js.common.HALF_PI;
1529            }
1530            else
1531            {  /* center of earth */
1532                Latitude = Proj4js.common.HALF_PI;
1533                Height = -this.b;
1534                return;
1535            }
1536        }
1537    }
1538    W2 = X*X + Y*Y;
1539    W = Math.sqrt(W2);
1540    T0 = Z * Proj4js.common.AD_C;
1541    S0 = Math.sqrt(T0 * T0 + W2);
1542    Sin_B0 = T0 / S0;
1543    Cos_B0 = W / S0;
1544    Sin3_B0 = Sin_B0 * Sin_B0 * Sin_B0;
1545    T1 = Z + this.b * this.ep2 * Sin3_B0;
1546    Sum = W - this.a * this.es * Cos_B0 * Cos_B0 * Cos_B0;
1547    S1 = Math.sqrt(T1*T1 + Sum * Sum);
1548    Sin_p1 = T1 / S1;
1549    Cos_p1 = Sum / S1;
1550    Rn = this.a / Math.sqrt(1.0 - this.es * Sin_p1 * Sin_p1);
1551    if (Cos_p1 >= Proj4js.common.COS_67P5)
1552    {
1553        Height = W / Cos_p1 - Rn;
1554    }
1555    else if (Cos_p1 <= -Proj4js.common.COS_67P5)
1556    {
1557        Height = W / -Cos_p1 - Rn;
1558    }
1559    else
1560    {
1561        Height = Z / Sin_p1 + Rn * (this.es - 1.0);
1562    }
1563    if (At_Pole == false)
1564    {
1565        Latitude = Math.atan(Sin_p1 / Cos_p1);
1566    }
1567
1568    p.x = Longitude;
1569    p.y = Latitude;
1570    p.z = Height;
1571    return p;
1572  }, // geocentric_to_geodetic_noniter()
1573
1574  /****************************************************************/
1575  // pj_geocentic_to_wgs84( p )
1576  //  p = point to transform in geocentric coordinates (x,y,z)
1577  geocentric_to_wgs84 : function ( p ) {
1578
1579    if( this.datum_type == Proj4js.common.PJD_3PARAM )
1580    {
1581      // if( x[io] == HUGE_VAL )
1582      //    continue;
1583      p.x += this.datum_params[0];
1584      p.y += this.datum_params[1];
1585      p.z += this.datum_params[2];
1586
1587    }
1588    else if (this.datum_type == Proj4js.common.PJD_7PARAM)
1589    {
1590      var Dx_BF =this.datum_params[0];
1591      var Dy_BF =this.datum_params[1];
1592      var Dz_BF =this.datum_params[2];
1593      var Rx_BF =this.datum_params[3];
1594      var Ry_BF =this.datum_params[4];
1595      var Rz_BF =this.datum_params[5];
1596      var M_BF  =this.datum_params[6];
1597      // if( x[io] == HUGE_VAL )
1598      //    continue;
1599      var x_out = M_BF*(       p.x - Rz_BF*p.y + Ry_BF*p.z) + Dx_BF;
1600      var y_out = M_BF*( Rz_BF*p.x +       p.y - Rx_BF*p.z) + Dy_BF;
1601      var z_out = M_BF*(-Ry_BF*p.x + Rx_BF*p.y +       p.z) + Dz_BF;
1602      p.x = x_out;
1603      p.y = y_out;
1604      p.z = z_out;
1605    }
1606  }, // cs_geocentric_to_wgs84
1607
1608  /****************************************************************/
1609  // pj_geocentic_from_wgs84()
1610  //  coordinate system definition,
1611  //  point to transform in geocentric coordinates (x,y,z)
1612  geocentric_from_wgs84 : function( p ) {
1613
1614    if( this.datum_type == Proj4js.common.PJD_3PARAM )
1615    {
1616      //if( x[io] == HUGE_VAL )
1617      //    continue;
1618      p.x -= this.datum_params[0];
1619      p.y -= this.datum_params[1];
1620      p.z -= this.datum_params[2];
1621
1622    }
1623    else if (this.datum_type == Proj4js.common.PJD_7PARAM)
1624    {
1625      var Dx_BF =this.datum_params[0];
1626      var Dy_BF =this.datum_params[1];
1627      var Dz_BF =this.datum_params[2];
1628      var Rx_BF =this.datum_params[3];
1629      var Ry_BF =this.datum_params[4];
1630      var Rz_BF =this.datum_params[5];
1631      var M_BF  =this.datum_params[6];
1632      var x_tmp = (p.x - Dx_BF) / M_BF;
1633      var y_tmp = (p.y - Dy_BF) / M_BF;
1634      var z_tmp = (p.z - Dz_BF) / M_BF;
1635      //if( x[io] == HUGE_VAL )
1636      //    continue;
1637
1638      p.x =        x_tmp + Rz_BF*y_tmp - Ry_BF*z_tmp;
1639      p.y = -Rz_BF*x_tmp +       y_tmp + Rx_BF*z_tmp;
1640      p.z =  Ry_BF*x_tmp - Rx_BF*y_tmp +       z_tmp;
1641    } //cs_geocentric_from_wgs84()
1642  }
1643});
1644
1645/** point object, nothing fancy, just allows values to be
1646    passed back and forth by reference rather than by value.
1647    Other point classes may be used as long as they have
1648    x and y properties, which will get modified in the transform method.
1649*/
1650Proj4js.Point = Proj4js.Class({
1651
1652    /**
1653     * Constructor: Proj4js.Point
1654     *
1655     * Parameters:
1656     * - x {float} or {Array} either the first coordinates component or
1657     *     the full coordinates
1658     * - y {float} the second component
1659     * - z {float} the third component, optional.
1660     */
1661    initialize : function(x,y,z) {
1662      if (typeof x == 'object') {
1663        this.x = x[0];
1664        this.y = x[1];
1665        this.z = x[2] || 0.0;
1666      } else if (typeof x == 'string' && typeof y == 'undefined') {
1667        var coords = x.split(',');
1668        this.x = parseFloat(coords[0]);
1669        this.y = parseFloat(coords[1]);
1670        this.z = parseFloat(coords[2]) || 0.0;
1671      } else {
1672        this.x = x;
1673        this.y = y;
1674        this.z = z || 0.0;
1675      }
1676    },
1677
1678    /**
1679     * APIMethod: clone
1680     * Build a copy of a Proj4js.Point object.
1681     *
1682     * Return:
1683     * {Proj4js}.Point the cloned point.
1684     */
1685    clone : function() {
1686      return new Proj4js.Point(this.x, this.y, this.z);
1687    },
1688
1689    /**
1690     * APIMethod: toString
1691     * Return a readable string version of the point
1692     *
1693     * Return:
1694     * {String} String representation of Proj4js.Point object.
1695     *           (ex. <i>"x=5,y=42"</i>)
1696     */
1697    toString : function() {
1698        return ("x=" + this.x + ",y=" + this.y);
1699    },
1700
1701    /**
1702     * APIMethod: toShortString
1703     * Return a short string version of the point.
1704     *
1705     * Return:
1706     * {String} Shortened String representation of Proj4js.Point object.
1707     *         (ex. <i>"5, 42"</i>)
1708     */
1709    toShortString : function() {
1710        return (this.x + ", " + this.y);
1711    }
1712});
1713
1714Proj4js.PrimeMeridian = {
1715    "greenwich": 0.0,               //"0dE",
1716    "lisbon":     -9.131906111111,   //"9d07'54.862\"W",
1717    "paris":       2.337229166667,   //"2d20'14.025\"E",
1718    "bogota":    -74.080916666667,  //"74d04'51.3\"W",
1719    "madrid":     -3.687938888889,  //"3d41'16.58\"W",
1720    "rome":       12.452333333333,  //"12d27'8.4\"E",
1721    "bern":        7.439583333333,  //"7d26'22.5\"E",
1722    "jakarta":   106.807719444444,  //"106d48'27.79\"E",
1723    "ferro":     -17.666666666667,  //"17d40'W",
1724    "brussels":    4.367975,        //"4d22'4.71\"E",
1725    "stockholm":  18.058277777778,  //"18d3'29.8\"E",
1726    "athens":     23.7163375,       //"23d42'58.815\"E",
1727    "oslo":       10.722916666667   //"10d43'22.5\"E"
1728};
1729
1730Proj4js.Ellipsoid = {
1731  "MERIT": {a:6378137.0, rf:298.257, ellipseName:"MERIT 1983"},
1732  "SGS85": {a:6378136.0, rf:298.257, ellipseName:"Soviet Geodetic System 85"},
1733  "GRS80": {a:6378137.0, rf:298.257222101, ellipseName:"GRS 1980(IUGG, 1980)"},
1734  "IAU76": {a:6378140.0, rf:298.257, ellipseName:"IAU 1976"},
1735  "airy": {a:6377563.396, b:6356256.910, ellipseName:"Airy 1830"},
1736  "APL4.": {a:6378137, rf:298.25, ellipseName:"Appl. Physics. 1965"},
1737  "NWL9D": {a:6378145.0, rf:298.25, ellipseName:"Naval Weapons Lab., 1965"},
1738  "mod_airy": {a:6377340.189, b:6356034.446, ellipseName:"Modified Airy"},
1739  "andrae": {a:6377104.43, rf:300.0, ellipseName:"Andrae 1876 (Den., Iclnd.)"},
1740  "aust_SA": {a:6378160.0, rf:298.25, ellipseName:"Australian Natl & S. Amer. 1969"},
1741  "GRS67": {a:6378160.0, rf:298.2471674270, ellipseName:"GRS 67(IUGG 1967)"},
1742  "bessel": {a:6377397.155, rf:299.1528128, ellipseName:"Bessel 1841"},
1743  "bess_nam": {a:6377483.865, rf:299.1528128, ellipseName:"Bessel 1841 (Namibia)"},
1744  "clrk66": {a:6378206.4, b:6356583.8, ellipseName:"Clarke 1866"},
1745  "clrk80": {a:6378249.145, rf:293.4663, ellipseName:"Clarke 1880 mod."},
1746  "CPM": {a:6375738.7, rf:334.29, ellipseName:"Comm. des Poids et Mesures 1799"},
1747  "delmbr": {a:6376428.0, rf:311.5, ellipseName:"Delambre 1810 (Belgium)"},
1748  "engelis": {a:6378136.05, rf:298.2566, ellipseName:"Engelis 1985"},
1749  "evrst30": {a:6377276.345, rf:300.8017, ellipseName:"Everest 1830"},
1750  "evrst48": {a:6377304.063, rf:300.8017, ellipseName:"Everest 1948"},
1751  "evrst56": {a:6377301.243, rf:300.8017, ellipseName:"Everest 1956"},
1752  "evrst69": {a:6377295.664, rf:300.8017, ellipseName:"Everest 1969"},
1753  "evrstSS": {a:6377298.556, rf:300.8017, ellipseName:"Everest (Sabah & Sarawak)"},
1754  "fschr60": {a:6378166.0, rf:298.3, ellipseName:"Fischer (Mercury Datum) 1960"},
1755  "fschr60m": {a:6378155.0, rf:298.3, ellipseName:"Fischer 1960"},
1756  "fschr68": {a:6378150.0, rf:298.3, ellipseName:"Fischer 1968"},
1757  "helmert": {a:6378200.0, rf:298.3, ellipseName:"Helmert 1906"},
1758  "hough": {a:6378270.0, rf:297.0, ellipseName:"Hough"},
1759  "intl": {a:6378388.0, rf:297.0, ellipseName:"International 1909 (Hayford)"},
1760  "kaula": {a:6378163.0, rf:298.24, ellipseName:"Kaula 1961"},
1761  "lerch": {a:6378139.0, rf:298.257, ellipseName:"Lerch 1979"},
1762  "mprts": {a:6397300.0, rf:191.0, ellipseName:"Maupertius 1738"},
1763  "new_intl": {a:6378157.5, b:6356772.2, ellipseName:"New International 1967"},
1764  "plessis": {a:6376523.0, rf:6355863.0, ellipseName:"Plessis 1817 (France)"},
1765  "krass": {a:6378245.0, rf:298.3, ellipseName:"Krassovsky, 1942"},
1766  "SEasia": {a:6378155.0, b:6356773.3205, ellipseName:"Southeast Asia"},
1767  "walbeck": {a:6376896.0, b:6355834.8467, ellipseName:"Walbeck"},
1768  "WGS60": {a:6378165.0, rf:298.3, ellipseName:"WGS 60"},
1769  "WGS66": {a:6378145.0, rf:298.25, ellipseName:"WGS 66"},
1770  "WGS72": {a:6378135.0, rf:298.26, ellipseName:"WGS 72"},
1771  "WGS84": {a:6378137.0, rf:298.257223563, ellipseName:"WGS 84"},
1772  "sphere": {a:6370997.0, b:6370997.0, ellipseName:"Normal Sphere (r=6370997)"}
1773};
1774
1775Proj4js.Datum = {
1776  "WGS84": {towgs84: "0,0,0", ellipse: "WGS84", datumName: "WGS84"},
1777  "GGRS87": {towgs84: "-199.87,74.79,246.62", ellipse: "GRS80", datumName: "Greek_Geodetic_Reference_System_1987"},
1778  "NAD83": {towgs84: "0,0,0", ellipse: "GRS80", datumName: "North_American_Datum_1983"},
1779  "NAD27": {nadgrids: "@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat", ellipse: "clrk66", datumName: "North_American_Datum_1927"},
1780  "potsdam": {towgs84: "606.0,23.0,413.0", ellipse: "bessel", datumName: "Potsdam Rauenberg 1950 DHDN"},
1781  "carthage": {towgs84: "-263.0,6.0,431.0", ellipse: "clark80", datumName: "Carthage 1934 Tunisia"},
1782  "hermannskogel": {towgs84: "653.0,-212.0,449.0", ellipse: "bessel", datumName: "Hermannskogel"},
1783  "ire65": {towgs84: "482.530,-130.596,564.557,-1.042,-0.214,-0.631,8.15", ellipse: "mod_airy", datumName: "Ireland 1965"},
1784  "nzgd49": {towgs84: "59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993", ellipse: "intl", datumName: "New Zealand Geodetic Datum 1949"},
1785  "OSGB36": {towgs84: "446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894", ellipse: "airy", datumName: "Airy 1830"}
1786};
1787
1788Proj4js.WGS84 = new Proj4js.Proj('WGS84');
1789Proj4js.Datum['OSB36'] = Proj4js.Datum['OSGB36']; //as returned from spatialreference.org
1790
1791//lookup table to go from the projection name in WKT to the Proj4js projection name
1792//build this out as required
1793Proj4js.wktProjections = {
1794  "Lambert Tangential Conformal Conic Projection": "lcc",
1795  "Mercator": "merc",
1796  "Popular Visualisation Pseudo Mercator": "merc",
1797  "Transverse_Mercator": "tmerc",
1798  "Transverse Mercator": "tmerc",
1799  "Lambert Azimuthal Equal Area": "laea",
1800  "Universal Transverse Mercator System": "utm"
1801};
1802
1803
1804/* ======================================================================
1805    projCode/aea.js
1806   ====================================================================== */
1807
1808/*******************************************************************************
1809NAME                     ALBERS CONICAL EQUAL AREA
1810
1811PURPOSE:        Transforms input longitude and latitude to Easting and Northing
1812                for the Albers Conical Equal Area projection.  The longitude
1813                and latitude must be in radians.  The Easting and Northing
1814                values will be returned in meters.
1815
1816PROGRAMMER              DATE
1817----------              ----
1818T. Mittan,              Feb, 1992
1819
1820ALGORITHM REFERENCES
1821
18221.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
1823    Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
1824    State Government Printing Office, Washington D.C., 1987.
1825
18262.  Snyder, John P. and Voxland, Philip M., "An Album of Map Projections",
1827    U.S. Geological Survey Professional Paper 1453 , United State Government
1828    Printing Office, Washington D.C., 1989.
1829*******************************************************************************/
1830
1831
1832Proj4js.Proj.aea = {
1833  init : function() {
1834
1835    if (Math.abs(this.lat1 + this.lat2) < Proj4js.common.EPSLN) {
1836       Proj4js.reportError("aeaInitEqualLatitudes");
1837       return;
1838    }
1839    this.temp = this.b / this.a;
1840    this.es = 1.0 - Math.pow(this.temp,2);
1841    this.e3 = Math.sqrt(this.es);
1842
1843    this.sin_po=Math.sin(this.lat1);
1844    this.cos_po=Math.cos(this.lat1);
1845    this.t1=this.sin_po;
1846    this.con = this.sin_po;
1847    this.ms1 = Proj4js.common.msfnz(this.e3,this.sin_po,this.cos_po);
1848    this.qs1 = Proj4js.common.qsfnz(this.e3,this.sin_po,this.cos_po);
1849
1850    this.sin_po=Math.sin(this.lat2);
1851    this.cos_po=Math.cos(this.lat2);
1852    this.t2=this.sin_po;
1853    this.ms2 = Proj4js.common.msfnz(this.e3,this.sin_po,this.cos_po);
1854    this.qs2 = Proj4js.common.qsfnz(this.e3,this.sin_po,this.cos_po);
1855
1856    this.sin_po=Math.sin(this.lat0);
1857    this.cos_po=Math.cos(this.lat0);
1858    this.t3=this.sin_po;
1859    this.qs0 = Proj4js.common.qsfnz(this.e3,this.sin_po,this.cos_po);
1860
1861    if (Math.abs(this.lat1 - this.lat2) > Proj4js.common.EPSLN) {
1862      this.ns0 = (this.ms1 * this.ms1 - this.ms2 *this.ms2)/ (this.qs2 - this.qs1);
1863    } else {
1864      this.ns0 = this.con;
1865    }
1866    this.c = this.ms1 * this.ms1 + this.ns0 * this.qs1;
1867    this.rh = this.a * Math.sqrt(this.c - this.ns0 * this.qs0)/this.ns0;
1868  },
1869
1870/* Albers Conical Equal Area forward equations--mapping lat,long to x,y
1871  -------------------------------------------------------------------*/
1872  forward: function(p){
1873
1874    var lon=p.x;
1875    var lat=p.y;
1876
1877    this.sin_phi=Math.sin(lat);
1878    this.cos_phi=Math.cos(lat);
1879
1880    var qs = Proj4js.common.qsfnz(this.e3,this.sin_phi,this.cos_phi);
1881    var rh1 =this.a * Math.sqrt(this.c - this.ns0 * qs)/this.ns0;
1882    var theta = this.ns0 * Proj4js.common.adjust_lon(lon - this.long0); 
1883    var x = rh1 * Math.sin(theta) + this.x0;
1884    var y = this.rh - rh1 * Math.cos(theta) + this.y0;
1885
1886    p.x = x; 
1887    p.y = y;
1888    return p;
1889  },
1890
1891
1892  inverse: function(p) {
1893    var rh1,qs,con,theta,lon,lat;
1894
1895    p.x -= this.x0;
1896    p.y = this.rh - p.y + this.y0;
1897    if (this.ns0 >= 0) {
1898      rh1 = Math.sqrt(p.x *p.x + p.y * p.y);
1899      con = 1.0;
1900    } else {
1901      rh1 = -Math.sqrt(p.x * p.x + p.y *p.y);
1902      con = -1.0;
1903    }
1904    theta = 0.0;
1905    if (rh1 != 0.0) {
1906      theta = Math.atan2(con * p.x, con * p.y);
1907    }
1908    con = rh1 * this.ns0 / this.a;
1909    qs = (this.c - con * con) / this.ns0;
1910    if (this.e3 >= 1e-10) {
1911      con = 1 - .5 * (1.0 -this.es) * Math.log((1.0 - this.e3) / (1.0 + this.e3))/this.e3;
1912      if (Math.abs(Math.abs(con) - Math.abs(qs)) > .0000000001 ) {
1913          lat = this.phi1z(this.e3,qs);
1914      } else {
1915          if (qs >= 0) {
1916             lat = .5 * PI;
1917          } else {
1918             lat = -.5 * PI;
1919          }
1920      }
1921    } else {
1922      lat = this.phi1z(e3,qs);
1923    }
1924
1925    lon = Proj4js.common.adjust_lon(theta/this.ns0 + this.long0);
1926    p.x = lon;
1927    p.y = lat;
1928    return p;
1929  },
1930 
1931/* Function to compute phi1, the latitude for the inverse of the
1932   Albers Conical Equal-Area projection.
1933-------------------------------------------*/
1934  phi1z: function (eccent,qs) {
1935    var con, com, dphi;
1936    var phi = Proj4js.common.asinz(.5 * qs);
1937    if (eccent < Proj4js.common.EPSLN) return phi;
1938   
1939    var eccnts = eccent * eccent; 
1940    for (var i = 1; i <= 25; i++) {
1941        sinphi = Math.sin(phi);
1942        cosphi = Math.cos(phi);
1943        con = eccent * sinphi; 
1944        com = 1.0 - con * con;
1945        dphi = .5 * com * com / cosphi * (qs / (1.0 - eccnts) - sinphi / com + .5 / eccent * Math.log((1.0 - con) / (1.0 + con)));
1946        phi = phi + dphi;
1947        if (Math.abs(dphi) <= 1e-7) return phi;
1948    }
1949    Proj4js.reportError("aea:phi1z:Convergence error");
1950    return null;
1951  }
1952 
1953};
1954
1955
1956
1957/* ======================================================================
1958    projCode/sterea.js
1959   ====================================================================== */
1960
1961
1962Proj4js.Proj.sterea = {
1963  dependsOn : 'gauss',
1964
1965  init : function() {
1966    Proj4js.Proj['gauss'].init.apply(this);
1967    if (!this.rc) {
1968      Proj4js.reportError("sterea:init:E_ERROR_0");
1969      return;
1970    }
1971    this.sinc0 = Math.sin(this.phic0);
1972    this.cosc0 = Math.cos(this.phic0);
1973    this.R2 = 2.0 * this.rc;
1974    if (!this.title) this.title = "Oblique Stereographic Alternative";
1975  },
1976
1977  forward : function(p) {
1978    p.x = Proj4js.common.adjust_lon(p.x-this.long0); /* adjust del longitude */
1979    Proj4js.Proj['gauss'].forward.apply(this, [p]);
1980    sinc = Math.sin(p.y);
1981    cosc = Math.cos(p.y);
1982    cosl = Math.cos(p.x);
1983    k = this.k0 * this.R2 / (1.0 + this.sinc0 * sinc + this.cosc0 * cosc * cosl);
1984    p.x = k * cosc * Math.sin(p.x);
1985    p.y = k * (this.cosc0 * sinc - this.sinc0 * cosc * cosl);
1986    p.x = this.a * p.x + this.x0;
1987    p.y = this.a * p.y + this.y0;
1988    return p;
1989  },
1990
1991  inverse : function(p) {
1992    var lon,lat;
1993    p.x = (p.x - this.x0) / this.a; /* descale and de-offset */
1994    p.y = (p.y - this.y0) / this.a;
1995
1996    p.x /= this.k0;
1997    p.y /= this.k0;
1998    if ( (rho = Math.sqrt(p.x*p.x + p.y*p.y)) ) {
1999      c = 2.0 * Math.atan2(rho, this.R2);
2000      sinc = Math.sin(c);
2001      cosc = Math.cos(c);
2002      lat = Math.asin(cosc * this.sinc0 + p.y * sinc * this.cosc0 / rho);
2003      lon = Math.atan2(p.x * sinc, rho * this.cosc0 * cosc - p.y * this.sinc0 * sinc);
2004    } else {
2005      lat = this.phic0;
2006      lon = 0.;
2007    }
2008
2009    p.x = lon;
2010    p.y = lat;
2011    Proj4js.Proj['gauss'].inverse.apply(this,[p]);
2012    p.x = Proj4js.common.adjust_lon(p.x + this.long0); /* adjust longitude to CM */
2013    return p;
2014  }
2015};
2016
2017/* ======================================================================
2018    projCode/poly.js
2019   ====================================================================== */
2020
2021/* Function to compute, phi4, the latitude for the inverse of the
2022   Polyconic projection.
2023------------------------------------------------------------*/
2024function phi4z (eccent,e0,e1,e2,e3,a,b,c,phi) {
2025        var sinphi, sin2ph, tanph, ml, mlp, con1, con2, con3, dphi, i;
2026
2027        phi = a;
2028        for (i = 1; i <= 15; i++) {
2029                sinphi = Math.sin(phi);
2030                tanphi = Math.tan(phi);
2031                c = tanphi * Math.sqrt (1.0 - eccent * sinphi * sinphi);
2032                sin2ph = Math.sin (2.0 * phi);
2033                /*
2034                ml = e0 * *phi - e1 * sin2ph + e2 * sin (4.0 *  *phi);
2035                mlp = e0 - 2.0 * e1 * cos (2.0 *  *phi) + 4.0 * e2 *  cos (4.0 *  *phi);
2036                */
2037                ml = e0 * phi - e1 * sin2ph + e2 * Math.sin (4.0 *  phi) - e3 * Math.sin (6.0 * phi);
2038                mlp = e0 - 2.0 * e1 * Math.cos (2.0 *  phi) + 4.0 * e2 * Math.cos (4.0 *  phi) - 6.0 * e3 * Math.cos (6.0 *  phi);
2039                con1 = 2.0 * ml + c * (ml * ml + b) - 2.0 * a *  (c * ml + 1.0);
2040                con2 = eccent * sin2ph * (ml * ml + b - 2.0 * a * ml) / (2.0 *c);
2041                con3 = 2.0 * (a - ml) * (c * mlp - 2.0 / sin2ph) - 2.0 * mlp;
2042                dphi = con1 / (con2 + con3);
2043                phi += dphi;
2044                if (Math.abs(dphi) <= .0000000001 ) return(phi);   
2045        }
2046        Proj4js.reportError("phi4z: No convergence");
2047        return null;
2048}
2049
2050
2051/* Function to compute the constant e4 from the input of the eccentricity
2052   of the spheroid, x.  This constant is used in the Polar Stereographic
2053   projection.
2054--------------------------------------------------------------------*/
2055function e4fn(x) {
2056        var con, com;
2057        con = 1.0 + x;
2058        com = 1.0 - x;
2059        return (Math.sqrt((Math.pow(con,con))*(Math.pow(com,com))));
2060}
2061
2062
2063
2064
2065
2066/*******************************************************************************
2067NAME                             POLYCONIC
2068
2069PURPOSE:        Transforms input longitude and latitude to Easting and
2070                Northing for the Polyconic projection.  The
2071                longitude and latitude must be in radians.  The Easting
2072                and Northing values will be returned in meters.
2073
2074PROGRAMMER              DATE
2075----------              ----
2076T. Mittan               Mar, 1993
2077
2078ALGORITHM REFERENCES
2079
20801.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
2081    Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
2082    State Government Printing Office, Washington D.C., 1987.
2083
20842.  Snyder, John P. and Voxland, Philip M., "An Album of Map Projections",
2085    U.S. Geological Survey Professional Paper 1453 , United State Government
2086    Printing Office, Washington D.C., 1989.
2087*******************************************************************************/
2088
2089Proj4js.Proj.poly = {
2090
2091        /* Initialize the POLYCONIC projection
2092          ----------------------------------*/
2093        init: function() {
2094                var temp;                       /* temporary variable           */
2095                if (this.lat0=0) this.lat0=90;//this.lat0 ca
2096
2097                /* Place parameters in static storage for common use
2098                  -------------------------------------------------*/
2099                this.temp = this.b / this.a;
2100                this.es = 1.0 - Math.pow(this.temp,2);// devait etre dans tmerc.js mais n y est pas donc je commente sinon retour de valeurs nulles
2101                this.e = Math.sqrt(this.es);
2102                this.e0 = Proj4js.common.e0fn(this.es);
2103                this.e1 = Proj4js.common.e1fn(this.es);
2104                this.e2 = Proj4js.common.e2fn(this.es);
2105                this.e3 = Proj4js.common.e3fn(this.es);
2106                this.ml0 = Proj4js.common.mlfn(this.e0, this.e1,this.e2, this.e3, this.lat0);//si que des zeros le calcul ne se fait pas
2107                //if (!this.ml0) {this.ml0=0;}
2108        },
2109
2110
2111        /* Polyconic forward equations--mapping lat,long to x,y
2112          ---------------------------------------------------*/
2113        forward: function(p) {
2114                var sinphi, cosphi;     /* sin and cos value                            */
2115                var al;                         /* temporary values                             */
2116                var c;                          /* temporary values                             */
2117                var con, ml;            /* cone constant, small m                       */
2118                var ms;                         /* small m                                      */
2119                var x,y;
2120
2121                var lon=p.x;
2122                var lat=p.y;   
2123
2124                con = Proj4js.common.adjust_lon(lon - this.long0);
2125                if (Math.abs(lat) <= .0000001) {
2126                        x = this.x0 + this.a * con;
2127                        y = this.y0 - this.a * this.ml0;
2128                } else {
2129                        sinphi = Math.sin(lat);
2130                        cosphi = Math.cos(lat);   
2131
2132                        ml = Proj4js.common.mlfn(this.e0, this.e1, this.e2, this.e3, lat);
2133                        ms = Proj4js.common.msfnz(this.e,sinphi,cosphi);
2134                        con = sinphi;
2135                        x = this.x0 + this.a * ms * Math.sin(con)/sinphi;
2136                        y = this.y0 + this.a * (ml - this.ml0 + ms * (1.0 - Math.cos(con))/sinphi);
2137                }
2138
2139                p.x=x;
2140                p.y=y;   
2141                return p;
2142        },
2143
2144
2145        /* Inverse equations
2146        -----------------*/
2147        inverse: function(p) {
2148                var sin_phi, cos_phi;   /* sin and cos value                            */
2149                var al;                                 /* temporary values                             */
2150                var b;                                  /* temporary values                             */
2151                var c;                                  /* temporary values                             */
2152                var con, ml;                    /* cone constant, small m                       */
2153                var iflg;                               /* error flag                                   */
2154                var lon,lat;
2155                p.x -= this.x0;
2156                p.y -= this.y0;
2157                al = this.ml0 + p.y/this.a;
2158                iflg = 0;
2159
2160                if (Math.abs(al) <= .0000001) {
2161                        lon = p.x/this.a + this.long0;
2162                        lat = 0.0;
2163                } else {
2164                        b = al * al + (p.x/this.a) * (p.x/this.a);
2165                        iflg = phi4z(this.es,this.e0,this.e1,this.e2,this.e3,this.al,b,c,lat);
2166                        if (iflg != 1) return(iflg);
2167                        lon = Proj4js.common.adjust_lon((Proj4js.common.asinz(p.x * c / this.a) / Math.sin(lat)) + this.long0);
2168                }
2169
2170                p.x=lon;
2171                p.y=lat;
2172                return p;
2173        }
2174};
2175
2176
2177
2178/* ======================================================================
2179    projCode/equi.js
2180   ====================================================================== */
2181
2182/*******************************************************************************
2183NAME                             EQUIRECTANGULAR
2184
2185PURPOSE:        Transforms input longitude and latitude to Easting and
2186                Northing for the Equirectangular projection.  The
2187                longitude and latitude must be in radians.  The Easting
2188                and Northing values will be returned in meters.
2189
2190PROGRAMMER              DATE
2191----------              ----
2192T. Mittan               Mar, 1993
2193
2194ALGORITHM REFERENCES
2195
21961.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
2197    Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
2198    State Government Printing Office, Washington D.C., 1987.
2199
22002.  Snyder, John P. and Voxland, Philip M., "An Album of Map Projections",
2201    U.S. Geological Survey Professional Paper 1453 , United State Government
2202    Printing Office, Washington D.C., 1989.
2203*******************************************************************************/
2204Proj4js.Proj.equi = {
2205
2206  init: function() {
2207    if(!this.x0) this.x0=0;
2208    if(!this.y0) this.y0=0;
2209    if(!this.lat0) this.lat0=0;
2210    if(!this.long0) this.long0=0;
2211    ///this.t2;
2212  },
2213
2214
2215
2216/* Equirectangular forward equations--mapping lat,long to x,y
2217  ---------------------------------------------------------*/
2218  forward: function(p) {
2219
2220    var lon=p.x;                               
2221    var lat=p.y;                       
2222
2223    var dlon = Proj4js.common.adjust_lon(lon - this.long0);
2224    var x = this.x0 +this. a * dlon *Math.cos(this.lat0);
2225    var y = this.y0 + this.a * lat;
2226
2227    this.t1=x;
2228    this.t2=Math.cos(this.lat0);
2229    p.x=x;
2230    p.y=y;
2231    return p;
2232  },  //equiFwd()
2233
2234
2235
2236/* Equirectangular inverse equations--mapping x,y to lat/long
2237  ---------------------------------------------------------*/
2238  inverse: function(p) {
2239
2240    p.x -= this.x0;
2241    p.y -= this.y0;
2242    var lat = p.y /this. a;
2243
2244    if ( Math.abs(lat) > Proj4js.common.HALF_PI) {
2245        Proj4js.reportError("equi:Inv:DataError");
2246    }
2247    var lon = Proj4js.common.adjust_lon(this.long0 + p.x / (this.a * Math.cos(this.lat0)));
2248    p.x=lon;
2249    p.y=lat;
2250  }//equiInv()
2251};
2252
2253
2254/* ======================================================================
2255    projCode/merc.js
2256   ====================================================================== */
2257
2258/*******************************************************************************
2259NAME                            MERCATOR
2260
2261PURPOSE:        Transforms input longitude and latitude to Easting and
2262                Northing for the Mercator projection.  The
2263                longitude and latitude must be in radians.  The Easting
2264                and Northing values will be returned in meters.
2265
2266PROGRAMMER              DATE
2267----------              ----
2268D. Steinwand, EROS      Nov, 1991
2269T. Mittan               Mar, 1993
2270
2271ALGORITHM REFERENCES
2272
22731.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
2274    Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
2275    State Government Printing Office, Washington D.C., 1987.
2276
22772.  Snyder, John P. and Voxland, Philip M., "An Album of Map Projections",
2278    U.S. Geological Survey Professional Paper 1453 , United State Government
2279    Printing Office, Washington D.C., 1989.
2280*******************************************************************************/
2281
2282//static double r_major = a;               /* major axis                                */
2283//static double r_minor = b;               /* minor axis                                */
2284//static double lon_center = long0;        /* Center longitude (projection center) */
2285//static double lat_origin =  lat0;        /* center latitude                   */
2286//static double e,es;                      /* eccentricity constants            */
2287//static double m1;                            /* small value m                 */
2288//static double false_northing = y0;   /* y offset in meters                    */
2289//static double false_easting = x0;        /* x offset in meters                        */
2290//scale_fact = k0
2291
2292Proj4js.Proj.merc = {
2293  init : function() {
2294        //?this.temp = this.r_minor / this.r_major;
2295        //this.temp = this.b / this.a;
2296        //this.es = 1.0 - Math.sqrt(this.temp);
2297        //this.e = Math.sqrt( this.es );
2298        //?this.m1 = Math.cos(this.lat_origin) / (Math.sqrt( 1.0 - this.es * Math.sin(this.lat_origin) * Math.sin(this.lat_origin)));
2299        //this.m1 = Math.cos(0.0) / (Math.sqrt( 1.0 - this.es * Math.sin(0.0) * Math.sin(0.0)));
2300    if (this.lat_ts) {
2301      if (this.sphere) {
2302        this.k0 = Math.cos(this.lat_ts);
2303      } else {
2304        this.k0 = Proj4js.common.msfnz(this.es, Math.sin(this.lat_ts), Math.cos(this.lat_ts));
2305      }
2306    }
2307  },
2308
2309/* Mercator forward equations--mapping lat,long to x,y
2310  --------------------------------------------------*/
2311
2312  forward : function(p) {       
2313    //alert("ll2m coords : "+coords);
2314    var lon = p.x;
2315    var lat = p.y;
2316    // convert to radians
2317    if ( lat*Proj4js.common.R2D > 90.0 && 
2318          lat*Proj4js.common.R2D < -90.0 && 
2319          lon*Proj4js.common.R2D > 180.0 && 
2320          lon*Proj4js.common.R2D < -180.0) {
2321      Proj4js.reportError("merc:forward: llInputOutOfRange: "+ lon +" : " + lat);
2322      return null;
2323    }
2324
2325    var x,y;
2326    if(Math.abs( Math.abs(lat) - Proj4js.common.HALF_PI)  <= Proj4js.common.EPSLN) {
2327      Proj4js.reportError("merc:forward: ll2mAtPoles");
2328      return null;
2329    } else {
2330      if (this.sphere) {
2331        x = this.x0 + this.a * this.k0 * Proj4js.common.adjust_lon(lon - this.long0);
2332        y = this.y0 + this.a * this.k0 * Math.log(Math.tan(Proj4js.common.FORTPI + 0.5*lat));
2333      } else {
2334        var sinphi = Math.sin(lat);
2335        var ts = Proj4js.common.tsfnz(this.e,lat,sinphi);
2336        x = this.x0 + this.a * this.k0 * Proj4js.common.adjust_lon(lon - this.long0);
2337        y = this.y0 - this.a * this.k0 * Math.log(ts);
2338      }
2339      p.x = x; 
2340      p.y = y;
2341      return p;
2342    }
2343  },
2344
2345
2346  /* Mercator inverse equations--mapping x,y to lat/long
2347  --------------------------------------------------*/
2348  inverse : function(p) {       
2349
2350    var x = p.x - this.x0;
2351    var y = p.y - this.y0;
2352    var lon,lat;
2353
2354    if (this.sphere) {
2355      lat = Proj4js.common.HALF_PI - 2.0 * Math.atan(Math.exp(-y / this.a * this.k0));
2356    } else {
2357      var ts = Math.exp(-y / (this.a * this.k0));
2358      lat = Proj4js.common.phi2z(this.e,ts);
2359      if(lat == -9999) {
2360        Proj4js.reportError("merc:inverse: lat = -9999");
2361        return null;
2362      }
2363    }
2364    lon = Proj4js.common.adjust_lon(this.long0+ x / (this.a * this.k0));
2365
2366    p.x = lon;
2367    p.y = lat;
2368    return p;
2369  }
2370};
2371
2372
2373/* ======================================================================
2374    projCode/utm.js
2375   ====================================================================== */
2376
2377/*******************************************************************************
2378NAME                            TRANSVERSE MERCATOR
2379
2380PURPOSE:        Transforms input longitude and latitude to Easting and
2381                Northing for the Transverse Mercator projection.  The
2382                longitude and latitude must be in radians.  The Easting
2383                and Northing values will be returned in meters.
2384
2385ALGORITHM REFERENCES
2386
23871.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
2388    Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
2389    State Government Printing Office, Washington D.C., 1987.
2390
23912.  Snyder, John P. and Voxland, Philip M., "An Album of Map Projections",
2392    U.S. Geological Survey Professional Paper 1453 , United State Government
2393    Printing Office, Washington D.C., 1989.
2394*******************************************************************************/
2395
2396
2397/**
2398  Initialize Transverse Mercator projection
2399*/
2400
2401Proj4js.Proj.utm = {
2402  dependsOn : 'tmerc',
2403
2404  init : function() {
2405    if (!this.zone) {
2406      Proj4js.reportError("utm:init: zone must be specified for UTM");
2407      return;
2408    }
2409    this.lat0 = 0.0;
2410    this.long0 = ((6 * Math.abs(this.zone)) - 183) * Proj4js.common.D2R;
2411    this.x0 = 500000.0;
2412    this.y0 = this.utmSouth ? 10000000.0 : 0.0;
2413    this.k0 = 0.9996;
2414
2415    Proj4js.Proj['tmerc'].init.apply(this);
2416    this.forward = Proj4js.Proj['tmerc'].forward;
2417    this.inverse = Proj4js.Proj['tmerc'].inverse;
2418  }
2419};
2420/* ======================================================================
2421    projCode/eqdc.js
2422   ====================================================================== */
2423
2424/*******************************************************************************
2425NAME                            EQUIDISTANT CONIC
2426
2427PURPOSE:        Transforms input longitude and latitude to Easting and Northing
2428                for the Equidistant Conic projection.  The longitude and
2429                latitude must be in radians.  The Easting and Northing values
2430                will be returned in meters.
2431
2432PROGRAMMER              DATE
2433----------              ----
2434T. Mittan               Mar, 1993
2435
2436ALGORITHM REFERENCES
2437
24381.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
2439    Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
2440    State Government Printing Office, Washington D.C., 1987.
2441
24422.  Snyder, John P. and Voxland, Philip M., "An Album of Map Projections",
2443    U.S. Geological Survey Professional Paper 1453 , United State Government
2444    Printing Office, Washington D.C., 1989.
2445*******************************************************************************/
2446
2447/* Variables common to all subroutines in this code file
2448  -----------------------------------------------------*/
2449
2450Proj4js.Proj.eqdc = {
2451
2452/* Initialize the Equidistant Conic projection
2453  ------------------------------------------*/
2454  init: function() {
2455
2456    /* Place parameters in static storage for common use
2457      -------------------------------------------------*/
2458
2459    if(!this.mode) this.mode=0;//chosen default mode
2460    this.temp = this.b / this.a;
2461    this.es = 1.0 - Math.pow(this.temp,2);
2462    this.e = Math.sqrt(this.es);
2463    this.e0 = Proj4js.common.e0fn(this.es);
2464    this.e1 = Proj4js.common.e1fn(this.es);
2465    this.e2 = Proj4js.common.e2fn(this.es);
2466    this.e3 = Proj4js.common.e3fn(this.es);
2467
2468    this.sinphi=Math.sin(this.lat1);
2469    this.cosphi=Math.cos(this.lat1);
2470
2471    this.ms1 = Proj4js.common.msfnz(this.e,this.sinphi,this.cosphi);
2472    this.ml1 = Proj4js.common.mlfn(this.e0, this.e1, this.e2,this.e3, this.lat1);
2473
2474    /* format B
2475    ---------*/
2476    if (this.mode != 0) {
2477      if (Math.abs(this.lat1 + this.lat2) < Proj4js.common.EPSLN) {
2478            Proj4js.reportError("eqdc:Init:EqualLatitudes");
2479            //return(81);
2480       }
2481       this.sinphi=Math.sin(this.lat2);
2482       this.cosphi=Math.cos(this.lat2);   
2483
2484       this.ms2 = Proj4js.common.msfnz(this.e,this.sinphi,this.cosphi);
2485       this.ml2 = Proj4js.common.mlfn(this.e0, this.e1, this.e2, this.e3, this.lat2);
2486       if (Math.abs(this.lat1 - this.lat2) >= Proj4js.common.EPSLN) {
2487         this.ns = (this.ms1 - this.ms2) / (this.ml2 - this.ml1);
2488       } else {
2489          this.ns = this.sinphi;
2490       }
2491    } else {
2492      this.ns = this.sinphi;
2493    }
2494    this.g = this.ml1 + this.ms1/this.ns;
2495    this.ml0 = Proj4js.common.mlfn(this.e0, this.e1,this. e2, this.e3, this.lat0);
2496    this.rh = this.a * (this.g - this.ml0);
2497  },
2498
2499
2500/* Equidistant Conic forward equations--mapping lat,long to x,y
2501  -----------------------------------------------------------*/
2502  forward: function(p) {
2503    var lon=p.x;
2504    var lat=p.y;
2505
2506    /* Forward equations
2507      -----------------*/
2508    var ml = Proj4js.common.mlfn(this.e0, this.e1, this.e2, this.e3, lat);
2509    var rh1 = this.a * (this.g - ml);
2510    var theta = this.ns * Proj4js.common.adjust_lon(lon - this.long0);
2511
2512    var x = this.x0  + rh1 * Math.sin(theta);
2513    var y = this.y0 + this.rh - rh1 * Math.cos(theta);
2514    p.x=x;
2515    p.y=y;
2516    return p;
2517  },
2518
2519/* Inverse equations
2520  -----------------*/
2521  inverse: function(p) {
2522    p.x -= this.x0;
2523    p.y  = this.rh - p.y + this.y0;
2524    var con, rh1;
2525    if (this.ns >= 0) {
2526       var rh1 = Math.sqrt(p.x *p.x + p.y * p.y); 
2527       var con = 1.0;
2528    } else {
2529       rh1 = -Math.sqrt(p.x *p. x +p. y * p.y); 
2530       con = -1.0;
2531    }
2532    var theta = 0.0;
2533    if (rh1 != 0.0) theta = Math.atan2(con *p.x, con *p.y);
2534    var ml = this.g - rh1 /this.a;
2535    var lat = this.phi3z(ml,this.e0,this.e1,this.e2,this.e3);
2536    var lon = Proj4js.common.adjust_lon(this.long0 + theta / this.ns);
2537
2538     p.x=lon;
2539     p.y=lat; 
2540     return p;
2541    },
2542   
2543/* Function to compute latitude, phi3, for the inverse of the Equidistant
2544   Conic projection.
2545-----------------------------------------------------------------*/
2546  phi3z: function(ml,e0,e1,e2,e3) {
2547    var phi;
2548    var dphi;
2549
2550    phi = ml;
2551    for (var i = 0; i < 15; i++) {
2552      dphi = (ml + e1 * Math.sin(2.0 * phi) - e2 * Math.sin(4.0 * phi) + e3 * Math.sin(6.0 * phi))/ e0 - phi;
2553      phi += dphi;
2554      if (Math.abs(dphi) <= .0000000001) {
2555        return phi;
2556      }
2557    }
2558    Proj4js.reportError("PHI3Z-CONV:Latitude failed to converge after 15 iterations");
2559    return null;
2560  }
2561
2562   
2563};
2564/* ======================================================================
2565    projCode/tmerc.js
2566   ====================================================================== */
2567
2568/*******************************************************************************
2569NAME                            TRANSVERSE MERCATOR
2570
2571PURPOSE:        Transforms input longitude and latitude to Easting and
2572                Northing for the Transverse Mercator projection.  The
2573                longitude and latitude must be in radians.  The Easting
2574                and Northing values will be returned in meters.
2575
2576ALGORITHM REFERENCES
2577
25781.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
2579    Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
2580    State Government Printing Office, Washington D.C., 1987.
2581
25822.  Snyder, John P. and Voxland, Philip M., "An Album of Map Projections",
2583    U.S. Geological Survey Professional Paper 1453 , United State Government
2584    Printing Office, Washington D.C., 1989.
2585*******************************************************************************/
2586
2587
2588/**
2589  Initialize Transverse Mercator projection
2590*/
2591
2592Proj4js.Proj.tmerc = {
2593  init : function() {
2594    this.e0 = Proj4js.common.e0fn(this.es);
2595    this.e1 = Proj4js.common.e1fn(this.es);
2596    this.e2 = Proj4js.common.e2fn(this.es);
2597    this.e3 = Proj4js.common.e3fn(this.es);
2598    this.ml0 = this.a * Proj4js.common.mlfn(this.e0, this.e1, this.e2, this.e3, this.lat0);
2599  },
2600
2601  /**
2602    Transverse Mercator Forward  - long/lat to x/y
2603    long/lat in radians
2604  */
2605  forward : function(p) {
2606    var lon = p.x;
2607    var lat = p.y;
2608
2609    var delta_lon = Proj4js.common.adjust_lon(lon - this.long0); // Delta longitude
2610    var con;    // cone constant
2611    var x, y;
2612    var sin_phi=Math.sin(lat);
2613    var cos_phi=Math.cos(lat);
2614
2615    if (this.sphere) {  /* spherical form */
2616      var b = cos_phi * Math.sin(delta_lon);
2617      if ((Math.abs(Math.abs(b) - 1.0)) < .0000000001)  {
2618        Proj4js.reportError("tmerc:forward: Point projects into infinity");
2619        return(93);
2620      } else {
2621        x = .5 * this.a * this.k0 * Math.log((1.0 + b)/(1.0 - b));
2622        con = Math.acos(cos_phi * Math.cos(delta_lon)/Math.sqrt(1.0 - b*b));
2623        if (lat < 0) con = - con;
2624        y = this.a * this.k0 * (con - this.lat0);
2625      }
2626    } else {
2627      var al  = cos_phi * delta_lon;
2628      var als = Math.pow(al,2);
2629      var c   = this.ep2 * Math.pow(cos_phi,2);
2630      var tq  = Math.tan(lat);
2631      var t   = Math.pow(tq,2);
2632      con = 1.0 - this.es * Math.pow(sin_phi,2);
2633      var n   = this.a / Math.sqrt(con);
2634      var ml  = this.a * Proj4js.common.mlfn(this.e0, this.e1, this.e2, this.e3, lat);
2635
2636      x = this.k0 * n * al * (1.0 + als / 6.0 * (1.0 - t + c + als / 20.0 * (5.0 - 18.0 * t + Math.pow(t,2) + 72.0 * c - 58.0 * this.ep2))) + this.x0;
2637      y = this.k0 * (ml - this.ml0 + n * tq * (als * (0.5 + als / 24.0 * (5.0 - t + 9.0 * c + 4.0 * Math.pow(c,2) + als / 30.0 * (61.0 - 58.0 * t + Math.pow(t,2) + 600.0 * c - 330.0 * this.ep2))))) + this.y0;
2638
2639    }
2640    p.x = x; p.y = y;
2641    return p;
2642  }, // tmercFwd()
2643
2644  /**
2645    Transverse Mercator Inverse  -  x/y to long/lat
2646  */
2647  inverse : function(p) {
2648    var con, phi;  /* temporary angles       */
2649    var delta_phi; /* difference between longitudes    */
2650    var i;
2651    var max_iter = 6;      /* maximun number of iterations */
2652    var lat, lon;
2653
2654    if (this.sphere) {   /* spherical form */
2655      var f = Math.exp(p.x/(this.a * this.k0));
2656      var g = .5 * (f - 1/f);
2657      var temp = this.lat0 + p.y/(this.a * this.k0);
2658      var h = Math.cos(temp);
2659      con = Math.sqrt((1.0 - h * h)/(1.0 + g * g));
2660      lat = Proj4js.common.asinz(con);
2661      if (temp < 0)
2662        lat = -lat;
2663      if ((g == 0) && (h == 0)) {
2664        lon = this.long0;
2665      } else {
2666        lon = Proj4js.common.adjust_lon(Math.atan2(g,h) + this.long0);
2667      }
2668    } else {    // ellipsoidal form
2669      var x = p.x - this.x0;
2670      var y = p.y - this.y0;
2671
2672      con = (this.ml0 + y / this.k0) / this.a;
2673      phi = con;
2674      for (i=0;true;i++) {
2675        delta_phi=((con + this.e1 * Math.sin(2.0*phi) - this.e2 * Math.sin(4.0*phi) + this.e3 * Math.sin(6.0*phi)) / this.e0) - phi;
2676        phi += delta_phi;
2677        if (Math.abs(delta_phi) <= Proj4js.common.EPSLN) break;
2678        if (i >= max_iter) {
2679          Proj4js.reportError("tmerc:inverse: Latitude failed to converge");
2680          return(95);
2681        }
2682      } // for()
2683      if (Math.abs(phi) < Proj4js.common.HALF_PI) {
2684        // sincos(phi, &sin_phi, &cos_phi);
2685        var sin_phi=Math.sin(phi);
2686        var cos_phi=Math.cos(phi);
2687        var tan_phi = Math.tan(phi);
2688        var c = this.ep2 * Math.pow(cos_phi,2);
2689        var cs = Math.pow(c,2);
2690        var t = Math.pow(tan_phi,2);
2691        var ts = Math.pow(t,2);
2692        con = 1.0 - this.es * Math.pow(sin_phi,2);
2693        var n = this.a / Math.sqrt(con);
2694        var r = n * (1.0 - this.es) / con;
2695        var d = x / (n * this.k0);
2696        var ds = Math.pow(d,2);
2697        lat = phi - (n * tan_phi * ds / r) * (0.5 - ds / 24.0 * (5.0 + 3.0 * t + 10.0 * c - 4.0 * cs - 9.0 * this.ep2 - ds / 30.0 * (61.0 + 90.0 * t + 298.0 * c + 45.0 * ts - 252.0 * this.ep2 - 3.0 * cs)));
2698        lon = Proj4js.common.adjust_lon(this.long0 + (d * (1.0 - ds / 6.0 * (1.0 + 2.0 * t + c - ds / 20.0 * (5.0 - 2.0 * c + 28.0 * t - 3.0 * cs + 8.0 * this.ep2 + 24.0 * ts))) / cos_phi));
2699      } else {
2700        lat = Proj4js.common.HALF_PI * Proj4js.common.sign(y);
2701        lon = this.long0;
2702      }
2703    }
2704    p.x = lon;
2705    p.y = lat;
2706    return p;
2707  } // tmercInv()
2708};
2709/* ======================================================================
2710    defs/GOOGLE.js
2711   ====================================================================== */
2712
2713Proj4js.defs["GOOGLE"]="+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs";
2714Proj4js.defs["EPSG:900913"]=Proj4js.defs["GOOGLE"];
2715/* ======================================================================
2716    projCode/gstmerc.js
2717   ====================================================================== */
2718
2719Proj4js.Proj.gstmerc = {
2720  init : function() {
2721
2722    // array of:  a, b, lon0, lat0, k0, x0, y0
2723      var temp= this.b / this.a;
2724      this.e= Math.sqrt(1.0 - temp*temp);
2725      this.lc= this.long0;
2726      this.rs= Math.sqrt(1.0+this.e*this.e*Math.pow(Math.cos(this.lat0),4.0)/(1.0-this.e*this.e));
2727      var sinz= Math.sin(this.lat0);
2728      var pc= Math.asin(sinz/this.rs);
2729      var sinzpc= Math.sin(pc);
2730      this.cp= Proj4js.common.latiso(0.0,pc,sinzpc)-this.rs*Proj4js.common.latiso(this.e,this.lat0,sinz);
2731      this.n2= this.k0*this.a*Math.sqrt(1.0-this.e*this.e)/(1.0-this.e*this.e*sinz*sinz);
2732      this.xs= this.x0;
2733      this.ys= this.y0-this.n2*pc;
2734
2735      if (!this.title) this.title = "Gauss Schreiber transverse mercator";
2736    },
2737
2738
2739    // forward equations--mapping lat,long to x,y
2740    // -----------------------------------------------------------------
2741    forward : function(p) {
2742
2743      var lon= p.x;
2744      var lat= p.y;
2745
2746      var L= this.rs*(lon-this.lc);
2747      var Ls= this.cp+(this.rs*Proj4js.common.latiso(this.e,lat,Math.sin(lat)));
2748      var lat1= Math.asin(Math.sin(L)/Proj4js.common.cosh(Ls));
2749      var Ls1= Proj4js.common.latiso(0.0,lat1,Math.sin(lat1));
2750      p.x= this.xs+(this.n2*Ls1);
2751      p.y= this.ys+(this.n2*Math.atan(Proj4js.common.sinh(Ls)/Math.cos(L)));
2752      return p;
2753    },
2754
2755  // inverse equations--mapping x,y to lat/long
2756  // -----------------------------------------------------------------
2757  inverse : function(p) {
2758
2759    var x= p.x;
2760    var y= p.y;
2761
2762    var L= Math.atan(Proj4js.common.sinh((x-this.xs)/this.n2)/Math.cos((y-this.ys)/this.n2));
2763    var lat1= Math.asin(Math.sin((y-this.ys)/this.n2)/Proj4js.common.cosh((x-this.xs)/this.n2));
2764    var LC= Proj4js.common.latiso(0.0,lat1,Math.sin(lat1));
2765    p.x= this.lc+L/this.rs;
2766    p.y= Proj4js.common.invlatiso(this.e,(LC-this.cp)/this.rs);
2767    return p;
2768  }
2769
2770};
2771/* ======================================================================
2772    projCode/ortho.js
2773   ====================================================================== */
2774
2775/*******************************************************************************
2776NAME                             ORTHOGRAPHIC
2777
2778PURPOSE:        Transforms input longitude and latitude to Easting and
2779                Northing for the Orthographic projection.  The
2780                longitude and latitude must be in radians.  The Easting
2781                and Northing values will be returned in meters.
2782
2783PROGRAMMER              DATE
2784----------              ----
2785T. Mittan               Mar, 1993
2786
2787ALGORITHM REFERENCES
2788
27891.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
2790    Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
2791    State Government Printing Office, Washington D.C., 1987.
2792
27932.  Snyder, John P. and Voxland, Philip M., "An Album of Map Projections",
2794    U.S. Geological Survey Professional Paper 1453 , United State Government
2795    Printing Office, Washington D.C., 1989.
2796*******************************************************************************/
2797
2798Proj4js.Proj.ortho = {
2799
2800  /* Initialize the Orthographic projection
2801    -------------------------------------*/
2802  init: function(def) {
2803    //double temp;                      /* temporary variable           */
2804
2805    /* Place parameters in static storage for common use
2806      -------------------------------------------------*/;
2807    this.sin_p14=Math.sin(this.lat0);
2808    this.cos_p14=Math.cos(this.lat0);   
2809  },
2810
2811
2812  /* Orthographic forward equations--mapping lat,long to x,y
2813    ---------------------------------------------------*/
2814  forward: function(p) {
2815    var sinphi, cosphi; /* sin and cos value                            */
2816    var dlon;           /* delta longitude value                        */
2817    var coslon;         /* cos of longitude                             */
2818    var ksp;            /* scale factor                                 */
2819    var g;             
2820    var lon=p.x;
2821    var lat=p.y;       
2822    /* Forward equations
2823      -----------------*/
2824    dlon = Proj4js.common.adjust_lon(lon - this.long0);
2825
2826    sinphi=Math.sin(lat);
2827    cosphi=Math.cos(lat);       
2828
2829    coslon = Math.cos(dlon);
2830    g = this.sin_p14 * sinphi + this.cos_p14 * cosphi * coslon;
2831    ksp = 1.0;
2832    if ((g > 0) || (Math.abs(g) <= Proj4js.common.EPSLN)) {
2833      var x = this.a * ksp * cosphi * Math.sin(dlon);
2834      var y = this.y0 + this.a * ksp * (this.cos_p14 * sinphi - this.sin_p14 * cosphi * coslon);
2835    } else {
2836      Proj4js.reportError("orthoFwdPointError");
2837    }
2838    p.x=x;
2839    p.y=y;
2840    return p;
2841  },
2842
2843
2844  inverse: function(p) {
2845    var rh;             /* height above ellipsoid                       */
2846    var z;              /* angle                                        */
2847    var sinz,cosz;      /* sin of z and cos of z                        */
2848    var temp;
2849    var con;
2850    var lon , lat;
2851    /* Inverse equations
2852      -----------------*/
2853    p.x -= this.x0;
2854    p.y -= this.y0;
2855    rh = Math.sqrt(p.x * p.x + p.y * p.y);
2856    if (rh > this.a + .0000001) {
2857      Proj4js.reportError("orthoInvDataError");
2858    }
2859    z = Proj4js.common.asinz(rh / this.a);
2860
2861    sinz=Math.sin(z);
2862    cosz=Math.cos(z);
2863
2864    lon = this.long0;
2865    if (Math.abs(rh) <= Proj4js.common.EPSLN) {
2866      lat = this.lat0; 
2867    }
2868    lat = Proj4js.common.asinz(cosz * this.sin_p14 + (p.y * sinz * this.cos_p14)/rh);
2869    con = Math.abs(this.lat0) - Proj4js.common.HALF_PI;
2870    if (Math.abs(con) <= Proj4js.common.EPSLN) {
2871       if (this.lat0 >= 0) {
2872          lon = Proj4js.common.adjust_lon(this.long0 + Math.atan2(p.x, -p.y));
2873       } else {
2874          lon = Proj4js.common.adjust_lon(this.long0 -Math.atan2(-p.x, p.y));
2875       }
2876    }
2877    con = cosz - this.sin_p14 * Math.sin(lat);
2878    p.x=lon;
2879    p.y=lat;
2880    return p;
2881  }
2882};
2883
2884
2885/* ======================================================================
2886    projCode/somerc.js
2887   ====================================================================== */
2888
2889/*******************************************************************************
2890NAME                       SWISS OBLIQUE MERCATOR
2891
2892PURPOSE:        Swiss projection.
2893WARNING:  X and Y are inverted (weird) in the swiss coordinate system. Not
2894   here, since we want X to be horizontal and Y vertical.
2895
2896ALGORITHM REFERENCES
28971. "Formules et constantes pour le Calcul pour la
2898 projection cylindrique conforme à axe oblique et pour la transformation entre
2899 des systÚmes de référence".
2900 http://www.swisstopo.admin.ch/internet/swisstopo/fr/home/topics/survey/sys/refsys/switzerland.parsysrelated1.31216.downloadList.77004.DownloadFile.tmp/swissprojectionfr.pdf
2901
2902*******************************************************************************/
2903
2904Proj4js.Proj.somerc = {
2905
2906  init: function() {
2907    var phy0 = this.lat0;
2908    this.lambda0 = this.long0;
2909    var sinPhy0 = Math.sin(phy0);
2910    var semiMajorAxis = this.a;
2911    var invF = this.rf;
2912    var flattening = 1 / invF;
2913    var e2 = 2 * flattening - Math.pow(flattening, 2);
2914    var e = this.e = Math.sqrt(e2);
2915    this.R = this.k0 * semiMajorAxis * Math.sqrt(1 - e2) / (1 - e2 * Math.pow(sinPhy0, 2.0));
2916    this.alpha = Math.sqrt(1 + e2 / (1 - e2) * Math.pow(Math.cos(phy0), 4.0));
2917    this.b0 = Math.asin(sinPhy0 / this.alpha);
2918    this.K = Math.log(Math.tan(Math.PI / 4.0 + this.b0 / 2.0))
2919            - this.alpha
2920            * Math.log(Math.tan(Math.PI / 4.0 + phy0 / 2.0))
2921            + this.alpha
2922            * e / 2
2923            * Math.log((1 + e * sinPhy0)
2924            / (1 - e * sinPhy0));
2925  },
2926
2927
2928  forward: function(p) {
2929    var Sa1 = Math.log(Math.tan(Math.PI / 4.0 - p.y / 2.0));
2930    var Sa2 = this.e / 2.0
2931            * Math.log((1 + this.e * Math.sin(p.y))
2932            / (1 - this.e * Math.sin(p.y)));
2933    var S = -this.alpha * (Sa1 + Sa2) + this.K;
2934
2935        // spheric latitude
2936    var b = 2.0 * (Math.atan(Math.exp(S)) - Math.PI / 4.0);
2937
2938        // spheric longitude
2939    var I = this.alpha * (p.x - this.lambda0);
2940
2941        // psoeudo equatorial rotation
2942    var rotI = Math.atan(Math.sin(I)
2943            / (Math.sin(this.b0) * Math.tan(b) +
2944               Math.cos(this.b0) * Math.cos(I)));
2945
2946    var rotB = Math.asin(Math.cos(this.b0) * Math.sin(b) -
2947                         Math.sin(this.b0) * Math.cos(b) * Math.cos(I));
2948
2949    p.y = this.R / 2.0
2950            * Math.log((1 + Math.sin(rotB)) / (1 - Math.sin(rotB)))
2951            + this.y0;
2952    p.x = this.R * rotI + this.x0;
2953    return p;
2954  },
2955
2956  inverse: function(p) {
2957    var Y = p.x - this.x0;
2958    var X = p.y - this.y0;
2959
2960    var rotI = Y / this.R;
2961    var rotB = 2 * (Math.atan(Math.exp(X / this.R)) - Math.PI / 4.0);
2962
2963    var b = Math.asin(Math.cos(this.b0) * Math.sin(rotB)
2964            + Math.sin(this.b0) * Math.cos(rotB) * Math.cos(rotI));
2965    var I = Math.atan(Math.sin(rotI)
2966            / (Math.cos(this.b0) * Math.cos(rotI) - Math.sin(this.b0)
2967            * Math.tan(rotB)));
2968
2969    var lambda = this.lambda0 + I / this.alpha;
2970
2971    var S = 0.0;
2972    var phy = b;
2973    var prevPhy = -1000.0;
2974    var iteration = 0;
2975    while (Math.abs(phy - prevPhy) > 0.0000001)
2976    {
2977      if (++iteration > 20)
2978      {
2979        Proj4js.reportError("omercFwdInfinity");
2980        return;
2981      }
2982      //S = Math.log(Math.tan(Math.PI / 4.0 + phy / 2.0));
2983      S = 1.0
2984              / this.alpha
2985              * (Math.log(Math.tan(Math.PI / 4.0 + b / 2.0)) - this.K)
2986              + this.e
2987              * Math.log(Math.tan(Math.PI / 4.0
2988              + Math.asin(this.e * Math.sin(phy))
2989              / 2.0));
2990      prevPhy = phy;
2991      phy = 2.0 * Math.atan(Math.exp(S)) - Math.PI / 2.0;
2992    }
2993
2994    p.x = lambda;
2995    p.y = phy;
2996    return p;
2997  }
2998};
2999/* ======================================================================
3000    projCode/stere.js
3001   ====================================================================== */
3002
3003
3004// Initialize the Stereographic projection
3005
3006Proj4js.Proj.stere = {
3007  ssfn_: function(phit, sinphi, eccen) {
3008        sinphi *= eccen;
3009        return (Math.tan (.5 * (Proj4js.common.HALF_PI + phit)) * Math.pow((1. - sinphi) / (1. + sinphi), .5 * eccen));
3010  },
3011  TOL:  1.e-8,
3012  NITER:        8,
3013  CONV: 1.e-10,
3014  S_POLE:       0,
3015  N_POLE:       1,
3016  OBLIQ:        2,
3017  EQUIT:        3,
3018
3019  init : function() {
3020        this.phits = this.lat_ts ? this.lat_ts : Proj4js.common.HALF_PI;
3021    var t = Math.abs(this.lat0);
3022        if ((Math.abs(t) - Proj4js.common.HALF_PI) < Proj4js.common.EPSLN) {
3023                this.mode = this.lat0 < 0. ? this.S_POLE : this.N_POLE;
3024        } else {
3025                this.mode = t > Proj4js.common.EPSLN ? this.OBLIQ : this.EQUIT;
3026    }
3027        this.phits = Math.abs(this.phits);
3028        if (this.es) {
3029                var X;
3030
3031                switch (this.mode) {
3032                case this.N_POLE:
3033                case this.S_POLE:
3034                        if (Math.abs(this.phits - Proj4js.common.HALF_PI) < Proj4js.common.EPSLN) {
3035                                this.akm1 = 2. * this.k0 / Math.sqrt(Math.pow(1+this.e,1+this.e)*Math.pow(1-this.e,1-this.e));
3036                        } else {
3037          t = Math.sin(this.phits);
3038                                this.akm1 = Math.cos(this.phits) / Proj4js.common.tsfnz(this.e, this.phits, t);
3039                                t *= this.e;
3040                                this.akm1 /= Math.sqrt(1. - t * t);
3041                        }
3042                        break;
3043                case this.EQUIT:
3044                        this.akm1 = 2. * this.k0;
3045                        break;
3046                case this.OBLIQ:
3047                        t = Math.sin(this.lat0);
3048                        X = 2. * Math.atan(this.ssfn_(this.lat0, t, this.e)) - Proj4js.common.HALF_PI;
3049                        t *= this.e;
3050                        this.akm1 = 2. * this.k0 * Math.cos(this.lat0) / Math.sqrt(1. - t * t);
3051                        this.sinX1 = Math.sin(X);
3052                        this.cosX1 = Math.cos(X);
3053                        break;
3054                }
3055        } else {
3056                switch (this.mode) {
3057                case this.OBLIQ:
3058                        this.sinph0 = Math.sin(this.lat0);
3059                        this.cosph0 = Math.cos(this.lat0);
3060                case this.EQUIT:
3061                        this.akm1 = 2. * this.k0;
3062                        break;
3063                case this.S_POLE:
3064                case this.N_POLE:
3065                        this.akm1 = Math.abs(this.phits - Proj4js.common.HALF_PI) >= Proj4js.common.EPSLN ?
3066                           Math.cos(this.phits) / Math.tan(Proj4js.common.FORTPI - .5 * this.phits) :
3067                           2. * this.k0 ;
3068                        break;
3069                }
3070        }
3071  }, 
3072
3073// Stereographic forward equations--mapping lat,long to x,y
3074  forward: function(p) {
3075    var lon = p.x;
3076    lon = Proj4js.common.adjust_lon(lon - this.long0);
3077    var lat = p.y;
3078    var x, y;
3079   
3080    if (this.sphere) {
3081        var  sinphi, cosphi, coslam, sinlam;
3082
3083        sinphi = Math.sin(lat);
3084        cosphi = Math.cos(lat);
3085        coslam = Math.cos(lon);
3086        sinlam = Math.sin(lon);
3087        switch (this.mode) {
3088        case this.EQUIT:
3089                y = 1. + cosphi * coslam;
3090                if (y <= Proj4js.common.EPSLN) {
3091          F_ERROR;
3092        }
3093        y = this.akm1 / y;
3094                x = y * cosphi * sinlam;
3095        y *= sinphi;
3096                break;
3097        case this.OBLIQ:
3098                y = 1. + this.sinph0 * sinphi + this.cosph0 * cosphi * coslam;
3099                if (y <= Proj4js.common.EPSLN) {
3100          F_ERROR;
3101        }
3102        y = this.akm1 / y;
3103                x = y * cosphi * sinlam;
3104                y *= this.cosph0 * sinphi - this.sinph0 * cosphi * coslam;
3105                break;
3106        case this.N_POLE:
3107                coslam = -coslam;
3108                lat = -lat;
3109        //Note  no break here so it conitnues through S_POLE
3110        case this.S_POLE:
3111                if (Math.abs(lat - Proj4js.common.HALF_PI) < this.TOL) {
3112          F_ERROR;
3113        }
3114        y = this.akm1 * Math.tan(Proj4js.common.FORTPI + .5 * lat);
3115                x = sinlam * y;
3116                y *= coslam;
3117                break;
3118        }
3119    } else {
3120        coslam = Math.cos(lon);
3121        sinlam = Math.sin(lon);
3122        sinphi = Math.sin(lat);
3123        if (this.mode == this.OBLIQ || this.mode == this.EQUIT) {
3124        X = 2. * Math.atan(this.ssfn_(lat, sinphi, this.e));
3125                sinX = Math.sin(X - Proj4js.common.HALF_PI);
3126                cosX = Math.cos(X);
3127        }
3128        switch (this.mode) {
3129        case this.OBLIQ:
3130                A = this.akm1 / (this.cosX1 * (1. + this.sinX1 * sinX + this.cosX1 * cosX * coslam));
3131                y = A * (this.cosX1 * sinX - this.sinX1 * cosX * coslam);
3132                x = A * cosX;
3133                break;
3134        case this.EQUIT:
3135                A = 2. * this.akm1 / (1. + cosX * coslam);
3136                y = A * sinX;
3137                x = A * cosX;
3138                break;
3139        case this.S_POLE:
3140                lat = -lat;
3141                coslam = - coslam;
3142                sinphi = -sinphi;
3143        case this.N_POLE:
3144                x = this.akm1 * Proj4js.common.tsfnz(this.e, lat, sinphi);
3145                y = - x * coslam;
3146                break;
3147        }
3148        x = x * sinlam;
3149    }
3150    p.x = x*this.a + this.x0;
3151    p.y = y*this.a + this.y0;
3152    return p;
3153  },
3154
3155
3156//* Stereographic inverse equations--mapping x,y to lat/long
3157  inverse: function(p) {
3158    var x = (p.x - this.x0)/this.a;   /* descale and de-offset */
3159    var y = (p.y - this.y0)/this.a;
3160    var lon, lat;
3161
3162    var cosphi, sinphi, tp=0.0, phi_l=0.0, rho, halfe=0.0, pi2=0.0;
3163    var i;
3164
3165    if (this.sphere) {
3166        var  c, rh, sinc, cosc;
3167
3168      rh = Math.sqrt(x*x + y*y);
3169      c = 2. * Math.atan(rh / this.akm1);
3170        sinc = Math.sin(c);
3171        cosc = Math.cos(c);
3172        lon = 0.;
3173        switch (this.mode) {
3174        case this.EQUIT:
3175                if (Math.abs(rh) <= Proj4js.common.EPSLN) {
3176                        lat = 0.;
3177                } else {
3178                        lat = Math.asin(y * sinc / rh);
3179        }
3180                if (cosc != 0. || x != 0.) lon = Math.atan2(x * sinc, cosc * rh);
3181                break;
3182        case this.OBLIQ:
3183                if (Math.abs(rh) <= Proj4js.common.EPSLN) {
3184                        lat = this.phi0;
3185                } else {
3186                        lat = Math.asin(cosc * sinph0 + y * sinc * cosph0 / rh);
3187        }
3188        c = cosc - sinph0 * Math.sin(lat);
3189                if (c != 0. || x != 0.) {
3190                        lon = Math.atan2(x * sinc * cosph0, c * rh);
3191        }
3192                break;
3193        case this.N_POLE:
3194                y = -y;
3195        case this.S_POLE:
3196                if (Math.abs(rh) <= Proj4js.common.EPSLN) {
3197                        lat = this.phi0;
3198                } else {
3199                        lat = Math.asin(this.mode == this.S_POLE ? -cosc : cosc);
3200        }
3201                lon = (x == 0. && y == 0.) ? 0. : Math.atan2(x, y);
3202                break;
3203        }
3204        p.x = Proj4js.common.adjust_lon(lon + this.long0);
3205        p.y = lat;
3206    } else {
3207        rho = Math.sqrt(x*x + y*y);
3208        switch (this.mode) {
3209        case this.OBLIQ:
3210        case this.EQUIT:
3211        tp = 2. * Math.atan2(rho * this.cosX1 , this.akm1);
3212                cosphi = Math.cos(tp);
3213                sinphi = Math.sin(tp);
3214        if( rho == 0.0 ) {
3215                  phi_l = Math.asin(cosphi * this.sinX1);
3216        } else {
3217                  phi_l = Math.asin(cosphi * this.sinX1 + (y * sinphi * this.cosX1 / rho));
3218        }
3219
3220                tp = Math.tan(.5 * (Proj4js.common.HALF_PI + phi_l));
3221                x *= sinphi;
3222                y = rho * this.cosX1 * cosphi - y * this.sinX1* sinphi;
3223                pi2 = Proj4js.common.HALF_PI;
3224                halfe = .5 * this.e;
3225                break;
3226        case this.N_POLE:
3227                y = -y;
3228        case this.S_POLE:
3229        tp = - rho / this.akm1;
3230                phi_l = Proj4js.common.HALF_PI - 2. * Math.atan(tp);
3231                pi2 = -Proj4js.common.HALF_PI;
3232                halfe = -.5 * this.e;
3233                break;
3234        }
3235        for (i = this.NITER; i--; phi_l = lat) { //check this
3236                sinphi = this.e * Math.sin(phi_l);
3237                lat = 2. * Math.atan(tp * Math.pow((1.+sinphi)/(1.-sinphi), halfe)) - pi2;
3238                if (Math.abs(phi_l - lat) < this.CONV) {
3239                        if (this.mode == this.S_POLE) lat = -lat;
3240                        lon = (x == 0. && y == 0.) ? 0. : Math.atan2(x, y);
3241          p.x = Proj4js.common.adjust_lon(lon + this.long0);
3242          p.y = lat;
3243                        return p;
3244                }
3245        }
3246    }
3247  }
3248}; 
3249/* ======================================================================
3250    projCode/nzmg.js
3251   ====================================================================== */
3252
3253/*******************************************************************************
3254NAME                            NEW ZEALAND MAP GRID
3255
3256PURPOSE:        Transforms input longitude and latitude to Easting and
3257                Northing for the New Zealand Map Grid projection.  The
3258                longitude and latitude must be in radians.  The Easting
3259                and Northing values will be returned in meters.
3260
3261
3262ALGORITHM REFERENCES
3263
32641.  Department of Land and Survey Technical Circular 1973/32
3265      http://www.linz.govt.nz/docs/miscellaneous/nz-map-definition.pdf
3266
32672.  OSG Technical Report 4.1
3268      http://www.linz.govt.nz/docs/miscellaneous/nzmg.pdf
3269
3270
3271IMPLEMENTATION NOTES
3272
3273The two references use different symbols for the calculated values. This
3274implementation uses the variable names similar to the symbols in reference [1].
3275
3276The alogrithm uses different units for delta latitude and delta longitude.
3277The delta latitude is assumed to be in units of seconds of arc x 10^-5.
3278The delta longitude is the usual radians. Look out for these conversions.
3279
3280The algorithm is described using complex arithmetic. There were three
3281options:
3282   * find and use a Javascript library for complex arithmetic
3283   * write my own complex library
3284   * expand the complex arithmetic by hand to simple arithmetic
3285
3286This implementation has expanded the complex multiplication operations
3287into parallel simple arithmetic operations for the real and imaginary parts.
3288The imaginary part is way over to the right of the display; this probably
3289violates every coding standard in the world, but, to me, it makes it much
3290more obvious what is going on.
3291
3292The following complex operations are used:
3293   - addition
3294   - multiplication
3295   - division
3296   - complex number raised to integer power
3297   - summation
3298
3299A summary of complex arithmetic operations:
3300   (from http://en.wikipedia.org/wiki/Complex_arithmetic)
3301   addition:       (a + bi) + (c + di) = (a + c) + (b + d)i
3302   subtraction:    (a + bi) - (c + di) = (a - c) + (b - d)i
3303   multiplication: (a + bi) x (c + di) = (ac - bd) + (bc + ad)i
3304   division:       (a + bi) / (c + di) = [(ac + bd)/(cc + dd)] + [(bc - ad)/(cc + dd)]i
3305
3306The algorithm needs to calculate summations of simple and complex numbers. This is
3307implemented using a for-loop, pre-loading the summed value to zero.
3308
3309The algorithm needs to calculate theta^2, theta^3, etc while doing a summation.
3310There are three possible implementations:
3311   - use Math.pow in the summation loop - except for complex numbers
3312   - precalculate the values before running the loop
3313   - calculate theta^n = theta^(n-1) * theta during the loop
3314This implementation uses the third option for both real and complex arithmetic.
3315
3316For example
3317   psi_n = 1;
3318   sum = 0;
3319   for (n = 1; n <=6; n++) {
3320      psi_n1 = psi_n * psi;       // calculate psi^(n+1)
3321      psi_n = psi_n1;
3322      sum = sum + A[n] * psi_n;
3323   }
3324
3325
3326TEST VECTORS
3327
3328NZMG E, N:         2487100.638      6751049.719     metres
3329NZGD49 long, lat:      172.739194       -34.444066  degrees
3330
3331NZMG E, N:         2486533.395      6077263.661     metres
3332NZGD49 long, lat:      172.723106       -40.512409  degrees
3333
3334NZMG E, N:         2216746.425      5388508.765     metres
3335NZGD49 long, lat:      169.172062       -46.651295  degrees
3336
3337Note that these test vectors convert from NZMG metres to lat/long referenced
3338to NZGD49, not the more usual WGS84. The difference is about 70m N/S and about
333910m E/W.
3340
3341These test vectors are provided in reference [1]. Many more test
3342vectors are available in
3343   http://www.linz.govt.nz/docs/topography/topographicdata/placenamesdatabase/nznamesmar08.zip
3344which is a catalog of names on the 260-series maps.
3345
3346
3347EPSG CODES
3348
3349NZMG     EPSG:27200
3350NZGD49   EPSG:4272
3351
3352http://spatialreference.org/ defines these as
3353  Proj4js.defs["EPSG:4272"] = "+proj=longlat +ellps=intl +datum=nzgd49 +no_defs ";
3354  Proj4js.defs["EPSG:27200"] = "+proj=nzmg +lat_0=-41 +lon_0=173 +x_0=2510000 +y_0=6023150 +ellps=intl +datum=nzgd49 +units=m +no_defs ";
3355
3356
3357LICENSE
3358  Copyright: Stephen Irons 2008
3359  Released under terms of the LGPL as per: http://www.gnu.org/copyleft/lesser.html
3360
3361*******************************************************************************/
3362
3363
3364/**
3365  Initialize New Zealand Map Grip projection
3366*/
3367
3368Proj4js.Proj.nzmg = {
3369
3370  /**
3371   * iterations: Number of iterations to refine inverse transform.
3372   *     0 -> km accuracy
3373   *     1 -> m accuracy -- suitable for most mapping applications
3374   *     2 -> mm accuracy
3375   */
3376  iterations: 1,
3377
3378  init : function() {
3379    this.A = new Array();
3380    this.A[1]  = +0.6399175073;
3381    this.A[2]  = -0.1358797613;
3382    this.A[3]  = +0.063294409;
3383    this.A[4]  = -0.02526853;
3384    this.A[5]  = +0.0117879;
3385    this.A[6]  = -0.0055161;
3386    this.A[7]  = +0.0026906;
3387    this.A[8]  = -0.001333;
3388    this.A[9]  = +0.00067;
3389    this.A[10] = -0.00034;
3390
3391    this.B_re = new Array();        this.B_im = new Array();
3392    this.B_re[1] = +0.7557853228;   this.B_im[1] =  0.0;
3393    this.B_re[2] = +0.249204646;    this.B_im[2] = +0.003371507;
3394    this.B_re[3] = -0.001541739;    this.B_im[3] = +0.041058560;
3395    this.B_re[4] = -0.10162907;     this.B_im[4] = +0.01727609;
3396    this.B_re[5] = -0.26623489;     this.B_im[5] = -0.36249218;
3397    this.B_re[6] = -0.6870983;      this.B_im[6] = -1.1651967;
3398
3399    this.C_re = new Array();        this.C_im = new Array();
3400    this.C_re[1] = +1.3231270439;   this.C_im[1] =  0.0;
3401    this.C_re[2] = -0.577245789;    this.C_im[2] = -0.007809598;
3402    this.C_re[3] = +0.508307513;    this.C_im[3] = -0.112208952;
3403    this.C_re[4] = -0.15094762;     this.C_im[4] = +0.18200602;
3404    this.C_re[5] = +1.01418179;     this.C_im[5] = +1.64497696;
3405    this.C_re[6] = +1.9660549;      this.C_im[6] = +2.5127645;
3406
3407    this.D = new Array();
3408    this.D[1] = +1.5627014243;
3409    this.D[2] = +0.5185406398;
3410    this.D[3] = -0.03333098;
3411    this.D[4] = -0.1052906;
3412    this.D[5] = -0.0368594;
3413    this.D[6] = +0.007317;
3414    this.D[7] = +0.01220;
3415    this.D[8] = +0.00394;
3416    this.D[9] = -0.0013;
3417  },
3418
3419  /**
3420    New Zealand Map Grid Forward  - long/lat to x/y
3421    long/lat in radians
3422  */
3423  forward : function(p) {
3424    var lon = p.x;
3425    var lat = p.y;
3426
3427    var delta_lat = lat - this.lat0;
3428    var delta_lon = lon - this.long0;
3429
3430    // 1. Calculate d_phi and d_psi    ...                          // and d_lambda
3431    // For this algorithm, delta_latitude is in seconds of arc x 10-5, so we need to scale to those units. Longitude is radians.
3432    var d_phi = delta_lat / Proj4js.common.SEC_TO_RAD * 1E-5;       var d_lambda = delta_lon;
3433    var d_phi_n = 1;  // d_phi^0
3434
3435    var d_psi = 0;
3436    for (n = 1; n <= 10; n++) {
3437      d_phi_n = d_phi_n * d_phi;
3438      d_psi = d_psi + this.A[n] * d_phi_n;
3439    }
3440
3441    // 2. Calculate theta
3442    var th_re = d_psi;                                              var th_im = d_lambda;
3443
3444    // 3. Calculate z
3445    var th_n_re = 1;                                                var th_n_im = 0;  // theta^0
3446    var th_n_re1;                                                   var th_n_im1;
3447
3448    var z_re = 0;                                                   var z_im = 0;
3449    for (n = 1; n <= 6; n++) {
3450      th_n_re1 = th_n_re*th_re - th_n_im*th_im;                     th_n_im1 = th_n_im*th_re + th_n_re*th_im;
3451      th_n_re = th_n_re1;                                           th_n_im = th_n_im1;
3452      z_re = z_re + this.B_re[n]*th_n_re - this.B_im[n]*th_n_im;    z_im = z_im + this.B_im[n]*th_n_re + this.B_re[n]*th_n_im;
3453    }
3454
3455    // 4. Calculate easting and northing
3456    x = (z_im * this.a) + this.x0;
3457    y = (z_re * this.a) + this.y0;
3458
3459    p.x = x; p.y = y;
3460
3461    return p;
3462  },
3463
3464
3465  /**
3466    New Zealand Map Grid Inverse  -  x/y to long/lat
3467  */
3468  inverse : function(p) {
3469
3470    var x = p.x;
3471    var y = p.y;
3472
3473    var delta_x = x - this.x0;
3474    var delta_y = y - this.y0;
3475
3476    // 1. Calculate z
3477    var z_re = delta_y / this.a;                                              var z_im = delta_x / this.a;
3478
3479    // 2a. Calculate theta - first approximation gives km accuracy
3480    var z_n_re = 1;                                                           var z_n_im = 0;  // z^0
3481    var z_n_re1;                                                              var z_n_im1;
3482
3483    var th_re = 0;                                                            var th_im = 0;
3484    for (n = 1; n <= 6; n++) {
3485      z_n_re1 = z_n_re*z_re - z_n_im*z_im;                                    z_n_im1 = z_n_im*z_re + z_n_re*z_im;
3486      z_n_re = z_n_re1;                                                       z_n_im = z_n_im1;
3487      th_re = th_re + this.C_re[n]*z_n_re - this.C_im[n]*z_n_im;              th_im = th_im + this.C_im[n]*z_n_re + this.C_re[n]*z_n_im;
3488    }
3489
3490    // 2b. Iterate to refine the accuracy of the calculation
3491    //        0 iterations gives km accuracy
3492    //        1 iteration gives m accuracy -- good enough for most mapping applications
3493    //        2 iterations bives mm accuracy
3494    for (i = 0; i < this.iterations; i++) {
3495       var th_n_re = th_re;                                                      var th_n_im = th_im;
3496       var th_n_re1;                                                             var th_n_im1;
3497
3498       var num_re = z_re;                                                        var num_im = z_im;
3499       for (n = 2; n <= 6; n++) {
3500         th_n_re1 = th_n_re*th_re - th_n_im*th_im;                               th_n_im1 = th_n_im*th_re + th_n_re*th_im;
3501         th_n_re = th_n_re1;                                                     th_n_im = th_n_im1;
3502         num_re = num_re + (n-1)*(this.B_re[n]*th_n_re - this.B_im[n]*th_n_im);  num_im = num_im + (n-1)*(this.B_im[n]*th_n_re + this.B_re[n]*th_n_im);
3503       }
3504
3505       th_n_re = 1;                                                              th_n_im = 0;
3506       var den_re = this.B_re[1];                                                var den_im = this.B_im[1];
3507       for (n = 2; n <= 6; n++) {
3508         th_n_re1 = th_n_re*th_re - th_n_im*th_im;                               th_n_im1 = th_n_im*th_re + th_n_re*th_im;
3509         th_n_re = th_n_re1;                                                     th_n_im = th_n_im1;
3510         den_re = den_re + n * (this.B_re[n]*th_n_re - this.B_im[n]*th_n_im);    den_im = den_im + n * (this.B_im[n]*th_n_re + this.B_re[n]*th_n_im);
3511       }
3512
3513       // Complex division
3514       var den2 = den_re*den_re + den_im*den_im;
3515       th_re = (num_re*den_re + num_im*den_im) / den2;                           th_im = (num_im*den_re - num_re*den_im) / den2;
3516    }
3517
3518    // 3. Calculate d_phi              ...                                    // and d_lambda
3519    var d_psi = th_re;                                                        var d_lambda = th_im;
3520    var d_psi_n = 1;  // d_psi^0
3521
3522    var d_phi = 0;
3523    for (n = 1; n <= 9; n++) {
3524       d_psi_n = d_psi_n * d_psi;
3525       d_phi = d_phi + this.D[n] * d_psi_n;
3526    }
3527
3528    // 4. Calculate latitude and longitude
3529    // d_phi is calcuated in second of arc * 10^-5, so we need to scale back to radians. d_lambda is in radians.
3530    var lat = this.lat0 + (d_phi * Proj4js.common.SEC_TO_RAD * 1E5);
3531    var lon = this.long0 +  d_lambda;
3532
3533    p.x = lon;
3534    p.y = lat;
3535
3536    return p;
3537  }
3538};
3539/* ======================================================================
3540    projCode/mill.js
3541   ====================================================================== */
3542
3543/*******************************************************************************
3544NAME                    MILLER CYLINDRICAL
3545
3546PURPOSE:        Transforms input longitude and latitude to Easting and
3547                Northing for the Miller Cylindrical projection.  The
3548                longitude and latitude must be in radians.  The Easting
3549                and Northing values will be returned in meters.
3550
3551PROGRAMMER              DATE           
3552----------              ----           
3553T. Mittan               March, 1993
3554
3555This function was adapted from the Lambert Azimuthal Equal Area projection
3556code (FORTRAN) in the General Cartographic Transformation Package software
3557which is available from the U.S. Geological Survey National Mapping Division.
3558 
3559ALGORITHM REFERENCES
3560
35611.  "New Equal-Area Map Projections for Noncircular Regions", John P. Snyder,
3562    The American Cartographer, Vol 15, No. 4, October 1988, pp. 341-355.
3563
35642.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
3565    Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
3566    State Government Printing Office, Washington D.C., 1987.
3567
35683.  "Software Documentation for GCTP General Cartographic Transformation
3569    Package", U.S. Geological Survey National Mapping Division, May 1982.
3570*******************************************************************************/
3571
3572Proj4js.Proj.mill = {
3573
3574/* Initialize the Miller Cylindrical projection
3575  -------------------------------------------*/
3576  init: function() {
3577    //no-op
3578  },
3579
3580
3581  /* Miller Cylindrical forward equations--mapping lat,long to x,y
3582    ------------------------------------------------------------*/
3583  forward: function(p) {
3584    var lon=p.x;
3585    var lat=p.y;
3586    /* Forward equations
3587      -----------------*/
3588    var dlon = Proj4js.common.adjust_lon(lon -this.long0);
3589    var x = this.x0 + this.a * dlon;
3590    var y = this.y0 + this.a * Math.log(Math.tan((Proj4js.common.PI / 4.0) + (lat / 2.5))) * 1.25;
3591
3592    p.x=x;
3593    p.y=y;
3594    return p;
3595  },//millFwd()
3596
3597  /* Miller Cylindrical inverse equations--mapping x,y to lat/long
3598    ------------------------------------------------------------*/
3599  inverse: function(p) {
3600    p.x -= this.x0;
3601    p.y -= this.y0;
3602
3603    var lon = Proj4js.common.adjust_lon(this.long0 + p.x /this.a);
3604    var lat = 2.5 * (Math.atan(Math.exp(0.8*p.y/this.a)) - Proj4js.common.PI / 4.0);
3605
3606    p.x=lon;
3607    p.y=lat;
3608    return p;
3609  }//millInv()
3610};
3611/* ======================================================================
3612    projCode/gnom.js
3613   ====================================================================== */
3614
3615/*****************************************************************************
3616NAME                             GNOMONIC
3617
3618PURPOSE:        Transforms input longitude and latitude to Easting and
3619                Northing for the Gnomonic Projection.
3620                Implementation based on the existing sterea and ortho
3621                implementations.
3622
3623PROGRAMMER              DATE
3624----------              ----
3625Richard Marsden         November 2009
3626
3627ALGORITHM REFERENCES
3628
36291.  Snyder, John P., "Flattening the Earth - Two Thousand Years of Map
3630    Projections", University of Chicago Press 1993
3631
36322.  Wolfram Mathworld "Gnomonic Projection"
3633    http://mathworld.wolfram.com/GnomonicProjection.html
3634    Accessed: 12th November 2009
3635******************************************************************************/
3636
3637Proj4js.Proj.gnom = {
3638
3639  /* Initialize the Gnomonic projection
3640    -------------------------------------*/
3641  init: function(def) {
3642
3643    /* Place parameters in static storage for common use
3644      -------------------------------------------------*/
3645    this.sin_p14=Math.sin(this.lat0);
3646    this.cos_p14=Math.cos(this.lat0);
3647    // Approximation for projecting points to the horizon (infinity)
3648    this.infinity_dist = 1000 * this.a;
3649    this.rc = 1;
3650  },
3651
3652
3653  /* Gnomonic forward equations--mapping lat,long to x,y
3654    ---------------------------------------------------*/
3655  forward: function(p) {
3656    var sinphi, cosphi; /* sin and cos value                            */
3657    var dlon;           /* delta longitude value                        */
3658    var coslon;         /* cos of longitude                             */
3659    var ksp;            /* scale factor                                 */
3660    var g;             
3661    var lon=p.x;
3662    var lat=p.y;       
3663    /* Forward equations
3664      -----------------*/
3665    dlon = Proj4js.common.adjust_lon(lon - this.long0);
3666
3667    sinphi=Math.sin(lat);
3668    cosphi=Math.cos(lat);       
3669
3670    coslon = Math.cos(dlon);
3671    g = this.sin_p14 * sinphi + this.cos_p14 * cosphi * coslon;
3672    ksp = 1.0;
3673    if ((g > 0) || (Math.abs(g) <= Proj4js.common.EPSLN)) {
3674      x = this.x0 + this.a * ksp * cosphi * Math.sin(dlon) / g;
3675      y = this.y0 + this.a * ksp * (this.cos_p14 * sinphi - this.sin_p14 * cosphi * coslon) / g;
3676    } else {
3677      Proj4js.reportError("orthoFwdPointError");
3678
3679      // Point is in the opposing hemisphere and is unprojectable
3680      // We still need to return a reasonable point, so we project
3681      // to infinity, on a bearing
3682      // equivalent to the northern hemisphere equivalent
3683      // This is a reasonable approximation for short shapes and lines that
3684      // straddle the horizon.
3685
3686      x = this.x0 + this.infinity_dist * cosphi * Math.sin(dlon);
3687      y = this.y0 + this.infinity_dist * (this.cos_p14 * sinphi - this.sin_p14 * cosphi * coslon);
3688
3689    }
3690    p.x=x;
3691    p.y=y;
3692    return p;
3693  },
3694
3695
3696  inverse: function(p) {
3697    var rh;             /* Rho */
3698    var z;              /* angle */
3699    var sinc, cosc;
3700    var c;
3701    var lon , lat;
3702
3703    /* Inverse equations
3704      -----------------*/
3705    p.x = (p.x - this.x0) / this.a;
3706    p.y = (p.y - this.y0) / this.a;
3707
3708    p.x /= this.k0;
3709    p.y /= this.k0;
3710
3711    if ( (rh = Math.sqrt(p.x * p.x + p.y * p.y)) ) {
3712      c = Math.atan2(rh, this.rc);
3713      sinc = Math.sin(c);
3714      cosc = Math.cos(c);
3715
3716      lat = Proj4js.common.asinz(cosc*this.sin_p14 + (p.y*sinc*this.cos_p14) / rh);
3717      lon = Math.atan2(p.x*sinc, rh*this.cos_p14*cosc - p.y*this.sin_p14*sinc);
3718      lon = Proj4js.common.adjust_lon(this.long0+lon);
3719    } else {
3720      lat = this.phic0;
3721      lon = 0.0;
3722    }
3723 
3724    p.x=lon;
3725    p.y=lat;
3726    return p;
3727  }
3728};
3729
3730
3731/* ======================================================================
3732    projCode/sinu.js
3733   ====================================================================== */
3734
3735/*******************************************************************************
3736NAME                            SINUSOIDAL
3737
3738PURPOSE:        Transforms input longitude and latitude to Easting and
3739                Northing for the Sinusoidal projection.  The
3740                longitude and latitude must be in radians.  The Easting
3741                and Northing values will be returned in meters.
3742
3743PROGRAMMER              DATE           
3744----------              ----           
3745D. Steinwand, EROS      May, 1991     
3746
3747This function was adapted from the Sinusoidal projection code (FORTRAN) in the
3748General Cartographic Transformation Package software which is available from
3749the U.S. Geological Survey National Mapping Division.
3750 
3751ALGORITHM REFERENCES
3752
37531.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
3754    Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
3755    State Government Printing Office, Washington D.C., 1987.
3756
37572.  "Software Documentation for GCTP General Cartographic Transformation
3758    Package", U.S. Geological Survey National Mapping Division, May 1982.
3759*******************************************************************************/
3760
3761Proj4js.Proj.sinu = {
3762
3763        /* Initialize the Sinusoidal projection
3764          ------------------------------------*/
3765        init: function() {
3766                /* Place parameters in static storage for common use
3767                  -------------------------------------------------*/
3768                this.R = 6370997.0; //Radius of earth
3769        },
3770
3771        /* Sinusoidal forward equations--mapping lat,long to x,y
3772        -----------------------------------------------------*/
3773        forward: function(p) {
3774                var x,y,delta_lon;     
3775                var lon=p.x;
3776                var lat=p.y;   
3777                /* Forward equations
3778                -----------------*/
3779                delta_lon = Proj4js.common.adjust_lon(lon - this.long0);
3780                x = this.R * delta_lon * Math.cos(lat) + this.x0;
3781                y = this.R * lat + this.y0;
3782
3783                p.x=x;
3784                p.y=y; 
3785                return p;
3786        },
3787
3788        inverse: function(p) {
3789                var lat,temp,lon;       
3790
3791                /* Inverse equations
3792                  -----------------*/
3793                p.x -= this.x0;
3794                p.y -= this.y0;
3795                lat = p.y / this.R;
3796                if (Math.abs(lat) > Proj4js.common.HALF_PI) {
3797                    Proj4js.reportError("sinu:Inv:DataError");
3798                }
3799                temp = Math.abs(lat) - Proj4js.common.HALF_PI;
3800                if (Math.abs(temp) > Proj4js.common.EPSLN) {
3801                        temp = this.long0+ p.x / (this.R *Math.cos(lat));
3802                        lon = Proj4js.common.adjust_lon(temp);
3803                } else {
3804                        lon = this.long0;
3805                }
3806                 
3807                p.x=lon;
3808                p.y=lat;
3809                return p;
3810        }
3811};
3812
3813
3814/* ======================================================================
3815    projCode/vandg.js
3816   ====================================================================== */
3817
3818/*******************************************************************************
3819NAME                    VAN DER GRINTEN
3820
3821PURPOSE:        Transforms input Easting and Northing to longitude and
3822                latitude for the Van der Grinten projection.  The
3823                Easting and Northing must be in meters.  The longitude
3824                and latitude values will be returned in radians.
3825
3826PROGRAMMER              DATE           
3827----------              ----           
3828T. Mittan               March, 1993
3829
3830This function was adapted from the Van Der Grinten projection code
3831(FORTRAN) in the General Cartographic Transformation Package software
3832which is available from the U.S. Geological Survey National Mapping Division.
3833 
3834ALGORITHM REFERENCES
3835
38361.  "New Equal-Area Map Projections for Noncircular Regions", John P. Snyder,
3837    The American Cartographer, Vol 15, No. 4, October 1988, pp. 341-355.
3838
38392.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
3840    Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
3841    State Government Printing Office, Washington D.C., 1987.
3842
38433.  "Software Documentation for GCTP General Cartographic Transformation
3844    Package", U.S. Geological Survey National Mapping Division, May 1982.
3845*******************************************************************************/
3846
3847Proj4js.Proj.vandg = {
3848
3849/* Initialize the Van Der Grinten projection
3850  ----------------------------------------*/
3851        init: function() {
3852                this.R = 6370997.0; //Radius of earth
3853        },
3854
3855        forward: function(p) {
3856
3857                var lon=p.x;
3858                var lat=p.y;   
3859
3860                /* Forward equations
3861                -----------------*/
3862                var dlon = Proj4js.common.adjust_lon(lon - this.long0);
3863                var x,y;
3864
3865                if (Math.abs(lat) <= Proj4js.common.EPSLN) {
3866                        x = this.x0  + this.R * dlon;
3867                        y = this.y0;
3868                }
3869                var theta = Proj4js.common.asinz(2.0 * Math.abs(lat / Proj4js.common.PI));
3870                if ((Math.abs(dlon) <= Proj4js.common.EPSLN) || (Math.abs(Math.abs(lat) - Proj4js.common.HALF_PI) <= Proj4js.common.EPSLN)) {
3871                        x = this.x0;
3872                        if (lat >= 0) {
3873                                y = this.y0 + Proj4js.common.PI * this.R * Math.tan(.5 * theta);
3874                        } else {
3875                                y = this.y0 + Proj4js.common.PI * this.R * - Math.tan(.5 * theta);
3876                        }
3877                        //  return(OK);
3878                }
3879                var al = .5 * Math.abs((Proj4js.common.PI / dlon) - (dlon / Proj4js.common.PI));
3880                var asq = al * al;
3881                var sinth = Math.sin(theta);
3882                var costh = Math.cos(theta);
3883
3884                var g = costh / (sinth + costh - 1.0);
3885                var gsq = g * g;
3886                var m = g * (2.0 / sinth - 1.0);
3887                var msq = m * m;
3888                var con = Proj4js.common.PI * this.R * (al * (g - msq) + Math.sqrt(asq * (g - msq) * (g - msq) - (msq + asq) * (gsq - msq))) / (msq + asq);
3889                if (dlon < 0) {
3890                 con = -con;
3891                }
3892                x = this.x0 + con;
3893                con = Math.abs(con / (Proj4js.common.PI * this.R));
3894                if (lat >= 0) {
3895                 y = this.y0 + Proj4js.common.PI * this.R * Math.sqrt(1.0 - con * con - 2.0 * al * con);
3896                } else {
3897                 y = this.y0 - Proj4js.common.PI * this.R * Math.sqrt(1.0 - con * con - 2.0 * al * con);
3898                }
3899                p.x = x;
3900                p.y = y;
3901                return p;
3902        },
3903
3904/* Van Der Grinten inverse equations--mapping x,y to lat/long
3905  ---------------------------------------------------------*/
3906        inverse: function(p) {
3907                var dlon;
3908                var xx,yy,xys,c1,c2,c3;
3909                var al,asq;
3910                var a1;
3911                var m1;
3912                var con;
3913                var th1;
3914                var d;
3915
3916                /* inverse equations
3917                -----------------*/
3918                p.x -= this.x0;
3919                p.y -= this.y0;
3920                con = Proj4js.common.PI * this.R;
3921                xx = p.x / con;
3922                yy =p.y / con;
3923                xys = xx * xx + yy * yy;
3924                c1 = -Math.abs(yy) * (1.0 + xys);
3925                c2 = c1 - 2.0 * yy * yy + xx * xx;
3926                c3 = -2.0 * c1 + 1.0 + 2.0 * yy * yy + xys * xys;
3927                d = yy * yy / c3 + (2.0 * c2 * c2 * c2 / c3 / c3 / c3 - 9.0 * c1 * c2 / c3 /c3) / 27.0;
3928                a1 = (c1 - c2 * c2 / 3.0 / c3) / c3;
3929                m1 = 2.0 * Math.sqrt( -a1 / 3.0);
3930                con = ((3.0 * d) / a1) / m1;
3931                if (Math.abs(con) > 1.0) {
3932                        if (con >= 0.0) {
3933                                con = 1.0;
3934                        } else {
3935                                con = -1.0;
3936                        }
3937                }
3938                th1 = Math.acos(con) / 3.0;
3939                if (p.y >= 0) {
3940                        lat = (-m1 *Math.cos(th1 + Proj4js.common.PI / 3.0) - c2 / 3.0 / c3) * Proj4js.common.PI;
3941                } else {
3942                        lat = -(-m1 * Math.cos(th1 + Proj4js.common.PI / 3.0) - c2 / 3.0 / c3) * Proj4js.common.PI;
3943                }
3944
3945                if (Math.abs(xx) < Proj4js.common.EPSLN) {
3946                        lon = this.long0;
3947                }
3948                lon = Proj4js.common.adjust_lon(this.long0 + Proj4js.common.PI * (xys - 1.0 + Math.sqrt(1.0 + 2.0 * (xx * xx - yy * yy) + xys * xys)) / 2.0 / xx);
3949
3950                p.x=lon;
3951                p.y=lat;
3952                return p;
3953        }
3954};
3955/* ======================================================================
3956    projCode/cea.js
3957   ====================================================================== */
3958
3959/*******************************************************************************
3960NAME                    LAMBERT CYLINDRICAL EQUAL AREA
3961
3962PURPOSE:        Transforms input longitude and latitude to Easting and
3963                Northing for the Lambert Cylindrical Equal Area projection.
3964                This class of projection includes the Behrmann and
3965                Gall-Peters Projections.  The
3966                longitude and latitude must be in radians.  The Easting
3967                and Northing values will be returned in meters.
3968
3969PROGRAMMER              DATE           
3970----------              ----
3971R. Marsden              August 2009
3972Winwaed Software Tech LLC, http://www.winwaed.com
3973
3974This function was adapted from the Miller Cylindrical Projection in the Proj4JS
3975library.
3976
3977Note: This implementation assumes a Spherical Earth. The (commented) code
3978has been included for the ellipsoidal forward transform, but derivation of
3979the ellispoidal inverse transform is beyond me. Note that most of the
3980Proj4JS implementations do NOT currently support ellipsoidal figures.
3981Therefore this is not seen as a problem - especially this lack of support
3982is explicitly stated here.
3983 
3984ALGORITHM REFERENCES
3985
39861.  "Cartographic Projection Procedures for the UNIX Environment -
3987     A User's Manual" by Gerald I. Evenden, USGS Open File Report 90-284
3988    and Release 4 Interim Reports (2003)
3989
39902.  Snyder, John P., "Flattening the Earth - Two Thousand Years of Map
3991    Projections", Univ. Chicago Press, 1993
3992*******************************************************************************/
3993
3994Proj4js.Proj.cea = {
3995
3996/* Initialize the Cylindrical Equal Area projection
3997  -------------------------------------------*/
3998  init: function() {
3999    //no-op
4000  },
4001
4002
4003  /* Cylindrical Equal Area forward equations--mapping lat,long to x,y
4004    ------------------------------------------------------------*/
4005  forward: function(p) {
4006    var lon=p.x;
4007    var lat=p.y;
4008    /* Forward equations
4009      -----------------*/
4010    dlon = Proj4js.common.adjust_lon(lon -this.long0);
4011    var x = this.x0 + this.a * dlon * Math.cos(this.lat_ts);
4012    var y = this.y0 + this.a * Math.sin(lat) / Math.cos(this.lat_ts);
4013   /* Elliptical Forward Transform
4014      Not implemented due to a lack of a matchign inverse function
4015    {
4016      var Sin_Lat = Math.sin(lat);
4017      var Rn = this.a * (Math.sqrt(1.0e0 - this.es * Sin_Lat * Sin_Lat ));
4018      x = this.x0 + this.a * dlon * Math.cos(this.lat_ts);
4019      y = this.y0 + Rn * Math.sin(lat) / Math.cos(this.lat_ts);
4020    }
4021   */
4022
4023
4024    p.x=x;
4025    p.y=y;
4026    return p;
4027  },//ceaFwd()
4028
4029  /* Cylindrical Equal Area inverse equations--mapping x,y to lat/long
4030    ------------------------------------------------------------*/
4031  inverse: function(p) {
4032    p.x -= this.x0;
4033    p.y -= this.y0;
4034
4035    var lon = Proj4js.common.adjust_lon( this.long0 + (p.x / this.a) / Math.cos(this.lat_ts) );
4036
4037    var lat = Math.asin( (p.y/this.a) * Math.cos(this.lat_ts) );
4038
4039    p.x=lon;
4040    p.y=lat;
4041    return p;
4042  }//ceaInv()
4043};
4044/* ======================================================================
4045    projCode/eqc.js
4046   ====================================================================== */
4047
4048/* similar to equi.js FIXME proj4 uses eqc */
4049Proj4js.Proj.eqc = {
4050  init : function() {
4051
4052      if(!this.x0) this.x0=0;
4053      if(!this.y0) this.y0=0;
4054      if(!this.lat0) this.lat0=0;
4055      if(!this.long0) this.long0=0;
4056      if(!this.lat_ts) this.lat_ts=0;
4057      if (!this.title) this.title = "Equidistant Cylindrical (Plate Carre)";
4058
4059      this.rc= Math.cos(this.lat_ts);
4060    },
4061
4062
4063    // forward equations--mapping lat,long to x,y
4064    // -----------------------------------------------------------------
4065    forward : function(p) {
4066
4067      var lon= p.x;
4068      var lat= p.y;
4069
4070      var dlon = Proj4js.common.adjust_lon(lon - this.long0);
4071      var dlat = Proj4js.common.adjust_lat(lat - this.lat0 );
4072      p.x= this.x0 + (this.a*dlon*this.rc);
4073      p.y= this.y0 + (this.a*dlat        );
4074      return p;
4075    },
4076
4077  // inverse equations--mapping x,y to lat/long
4078  // -----------------------------------------------------------------
4079  inverse : function(p) {
4080
4081    var x= p.x;
4082    var y= p.y;
4083
4084    p.x= Proj4js.common.adjust_lon(this.long0 + ((x - this.x0)/(this.a*this.rc)));
4085    p.y= Proj4js.common.adjust_lat(this.lat0  + ((y - this.y0)/(this.a        )));
4086    return p;
4087  }
4088
4089};
4090/* ======================================================================
4091    projCode/cass.js
4092   ====================================================================== */
4093
4094/*******************************************************************************
4095NAME                            CASSINI
4096
4097PURPOSE:        Transforms input longitude and latitude to Easting and
4098                Northing for the Cassini projection.  The
4099                longitude and latitude must be in radians.  The Easting
4100                and Northing values will be returned in meters.
4101    Ported from PROJ.4.
4102
4103
4104ALGORITHM REFERENCES
4105
41061.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
4107    Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
4108    State Government Printing Office, Washington D.C., 1987.
4109
41102.  Snyder, John P. and Voxland, Philip M., "An Album of Map Projections",
4111    U.S. Geological Survey Professional Paper 1453 , United State Government
4112*******************************************************************************/
4113
4114
4115//Proj4js.defs["EPSG:28191"] = "+proj=cass +lat_0=31.73409694444445 +lon_0=35.21208055555556 +x_0=170251.555 +y_0=126867.909 +a=6378300.789 +b=6356566.435 +towgs84=-275.722,94.7824,340.894,-8.001,-4.42,-11.821,1 +units=m +no_defs";
4116
4117// Initialize the Cassini projection
4118// -----------------------------------------------------------------
4119
4120Proj4js.Proj.cass = {
4121  init : function() {
4122    if (!this.sphere) {
4123      this.en = this.pj_enfn(this.es)
4124      this.m0 = this.pj_mlfn(this.lat0, Math.sin(this.lat0), Math.cos(this.lat0), this.en);
4125    }
4126  },
4127
4128  C1:   .16666666666666666666,
4129  C2:   .00833333333333333333,
4130  C3:   .04166666666666666666,
4131  C4:   .33333333333333333333,
4132  C5:   .06666666666666666666,
4133
4134
4135/* Cassini forward equations--mapping lat,long to x,y
4136  -----------------------------------------------------------------------*/
4137  forward: function(p) {
4138
4139    /* Forward equations
4140      -----------------*/
4141    var x,y;
4142    var lam=p.x;
4143    var phi=p.y;
4144    lam = Proj4js.common.adjust_lon(lam - this.long0);
4145   
4146    if (this.sphere) {
4147      x = Math.asin(Math.cos(phi) * Math.sin(lam));
4148      y = Math.atan2(Math.tan(phi) , Math.cos(lam)) - this.phi0;
4149    } else {
4150        //ellipsoid
4151      this.n = Math.sin(phi);
4152      this.c = Math.cos(phi);
4153      y = this.pj_mlfn(phi, this.n, this.c, this.en);
4154      this.n = 1./Math.sqrt(1. - this.es * this.n * this.n);
4155      this.tn = Math.tan(phi); 
4156      this.t = this.tn * this.tn;
4157      this.a1 = lam * this.c;
4158      this.c *= this.es * this.c / (1 - this.es);
4159      this.a2 = this.a1 * this.a1;
4160      x = this.n * this.a1 * (1. - this.a2 * this.t * (this.C1 - (8. - this.t + 8. * this.c) * this.a2 * this.C2));
4161      y -= this.m0 - this.n * this.tn * this.a2 * (.5 + (5. - this.t + 6. * this.c) * this.a2 * this.C3);
4162    }
4163   
4164    p.x = this.a*x + this.x0;
4165    p.y = this.a*y + this.y0;
4166    return p;
4167  },//cassFwd()
4168
4169/* Inverse equations
4170  -----------------*/
4171  inverse: function(p) {
4172    p.x -= this.x0;
4173    p.y -= this.y0;
4174    var x = p.x/this.a;
4175    var y = p.y/this.a;
4176   
4177    if (this.sphere) {
4178      this.dd = y + this.lat0;
4179      phi = Math.asin(Math.sin(this.dd) * Math.cos(x));
4180      lam = Math.atan2(Math.tan(x), Math.cos(this.dd));
4181    } else {
4182      /* ellipsoid */
4183      ph1 = this.pj_inv_mlfn(this.m0 + y, this.es, this.en);
4184      this.tn = Math.tan(ph1); 
4185      this.t = this.tn * this.tn;
4186      this.n = Math.sin(ph1);
4187      this.r = 1. / (1. - this.es * this.n * this.n);
4188      this.n = Math.sqrt(this.r);
4189      this.r *= (1. - this.es) * this.n;
4190      this.dd = x / this.n;
4191      this.d2 = this.dd * this.dd;
4192      phi = ph1 - (this.n * this.tn / this.r) * this.d2 * (.5 - (1. + 3. * this.t) * this.d2 * this.C3);
4193      lam = this.dd * (1. + this.t * this.d2 * (-this.C4 + (1. + 3. * this.t) * this.d2 * this.C5)) / Math.cos(ph1);
4194    }
4195    p.x = Proj4js.common.adjust_lon(this.long0+lam);
4196    p.y = phi;
4197    return p;
4198  },//lamazInv()
4199
4200
4201  //code from the PROJ.4 pj_mlfn.c file;  this may be useful for other projections
4202  pj_enfn: function(es) {
4203    en = new Array();
4204    en[0] = this.C00 - es * (this.C02 + es * (this.C04 + es * (this.C06 + es * this.C08)));
4205    en[1] = es * (this.C22 - es * (this.C04 + es * (this.C06 + es * this.C08)));
4206    var t = es * es;
4207    en[2] = t * (this.C44 - es * (this.C46 + es * this.C48));
4208    t *= es;
4209    en[3] = t * (this.C66 - es * this.C68);
4210    en[4] = t * es * this.C88;
4211    return en;
4212  },
4213 
4214  pj_mlfn: function(phi, sphi, cphi, en) {
4215    cphi *= sphi;
4216    sphi *= sphi;
4217    return(en[0] * phi - cphi * (en[1] + sphi*(en[2]+ sphi*(en[3] + sphi*en[4]))));
4218  },
4219 
4220  pj_inv_mlfn: function(arg, es, en) {
4221    k = 1./(1.-es);
4222    phi = arg;
4223    for (i = Proj4js.common.MAX_ITER; i ; --i) { /* rarely goes over 2 iterations */
4224      s = Math.sin(phi);
4225      t = 1. - es * s * s;
4226      //t = this.pj_mlfn(phi, s, Math.cos(phi), en) - arg;
4227      //phi -= t * (t * Math.sqrt(t)) * k;
4228      t = (this.pj_mlfn(phi, s, Math.cos(phi), en) - arg) * (t * Math.sqrt(t)) * k;
4229      phi -= t;
4230      if (Math.abs(t) < Proj4js.common.EPSLN)
4231        return phi;
4232    }
4233    Proj4js.reportError("cass:pj_inv_mlfn: Convergence error");
4234    return phi;
4235  },
4236
4237/* meridinal distance for ellipsoid and inverse
4238**      8th degree - accurate to < 1e-5 meters when used in conjuction
4239**              with typical major axis values.
4240**      Inverse determines phi to EPS (1e-11) radians, about 1e-6 seconds.
4241*/
4242  C00: 1.0,
4243  C02: .25,
4244  C04: .046875,
4245  C06: .01953125,
4246  C08: .01068115234375,
4247  C22: .75,
4248  C44: .46875,
4249  C46: .01302083333333333333,
4250  C48: .00712076822916666666,
4251  C66: .36458333333333333333,
4252  C68: .00569661458333333333,
4253  C88: .3076171875
4254
4255}
4256/* ======================================================================
4257    projCode/gauss.js
4258   ====================================================================== */
4259
4260
4261Proj4js.Proj.gauss = {
4262
4263  init : function() {
4264    sphi = Math.sin(this.lat0);
4265    cphi = Math.cos(this.lat0); 
4266    cphi *= cphi;
4267    this.rc = Math.sqrt(1.0 - this.es) / (1.0 - this.es * sphi * sphi);
4268    this.C = Math.sqrt(1.0 + this.es * cphi * cphi / (1.0 - this.es));
4269    this.phic0 = Math.asin(sphi / this.C);
4270    this.ratexp = 0.5 * this.C * this.e;
4271    this.K = Math.tan(0.5 * this.phic0 + Proj4js.common.FORTPI) / (Math.pow(Math.tan(0.5*this.lat0 + Proj4js.common.FORTPI), this.C) * Proj4js.common.srat(this.e*sphi, this.ratexp));
4272  },
4273
4274  forward : function(p) {
4275    var lon = p.x;
4276    var lat = p.y;
4277
4278    p.y = 2.0 * Math.atan( this.K * Math.pow(Math.tan(0.5 * lat + Proj4js.common.FORTPI), this.C) * Proj4js.common.srat(this.e * Math.sin(lat), this.ratexp) ) - Proj4js.common.HALF_PI;
4279    p.x = this.C * lon;
4280    return p;
4281  },
4282
4283  inverse : function(p) {
4284    var DEL_TOL = 1e-14;
4285    var lon = p.x / this.C;
4286    var lat = p.y;
4287    num = Math.pow(Math.tan(0.5 * lat + Proj4js.common.FORTPI)/this.K, 1./this.C);
4288    for (var i = Proj4js.common.MAX_ITER; i>0; --i) {
4289      lat = 2.0 * Math.atan(num * Proj4js.common.srat(this.e * Math.sin(p.y), -0.5 * this.e)) - Proj4js.common.HALF_PI;
4290      if (Math.abs(lat - p.y) < DEL_TOL) break;
4291      p.y = lat;
4292    }   
4293    /* convergence failed */
4294    if (!i) {
4295      Proj4js.reportError("gauss:inverse:convergence failed");
4296      return null;
4297    }
4298    p.x = lon;
4299    p.y = lat;
4300    return p;
4301  }
4302};
4303
4304/* ======================================================================
4305    projCode/omerc.js
4306   ====================================================================== */
4307
4308/*******************************************************************************
4309NAME                       OBLIQUE MERCATOR (HOTINE)
4310
4311PURPOSE:        Transforms input longitude and latitude to Easting and
4312                Northing for the Oblique Mercator projection.  The
4313                longitude and latitude must be in radians.  The Easting
4314                and Northing values will be returned in meters.
4315
4316PROGRAMMER              DATE
4317----------              ----
4318T. Mittan               Mar, 1993
4319
4320ALGORITHM REFERENCES
4321
43221.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
4323    Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
4324    State Government Printing Office, Washington D.C., 1987.
4325
43262.  Snyder, John P. and Voxland, Philip M., "An Album of Map Projections",
4327    U.S. Geological Survey Professional Paper 1453 , United State Government
4328    Printing Office, Washington D.C., 1989.
4329*******************************************************************************/
4330
4331Proj4js.Proj.omerc = {
4332
4333  /* Initialize the Oblique Mercator  projection
4334    ------------------------------------------*/
4335  init: function() {
4336    if (!this.mode) this.mode=0;
4337    if (!this.lon1)   {this.lon1=0;this.mode=1;}
4338    if (!this.lon2)   this.lon2=0;
4339    if (!this.lat2)    this.lat2=0;
4340
4341    /* Place parameters in static storage for common use
4342      -------------------------------------------------*/
4343    var temp = this.b/ this.a;
4344    var es = 1.0 - Math.pow(temp,2);
4345    var e = Math.sqrt(es);
4346
4347    this.sin_p20=Math.sin(this.lat0);
4348    this.cos_p20=Math.cos(this.lat0);
4349
4350    this.con = 1.0 - this.es * this.sin_p20 * this.sin_p20;
4351    this.com = Math.sqrt(1.0 - es);
4352    this.bl = Math.sqrt(1.0 + this.es * Math.pow(this.cos_p20,4.0)/(1.0 - es));
4353    this.al = this.a * this.bl * this.k0 * this.com / this.con;
4354    if (Math.abs(this.lat0) < Proj4js.common.EPSLN) {
4355       this.ts = 1.0;
4356       this.d = 1.0;
4357       this.el = 1.0;
4358    } else {
4359       this.ts = Proj4js.common.tsfnz(this.e,this.lat0,this.sin_p20);
4360       this.con = Math.sqrt(this.con);
4361       this.d = this.bl * this.com / (this.cos_p20 * this.con);
4362       if ((this.d * this.d - 1.0) > 0.0) {
4363          if (this.lat0 >= 0.0) {
4364             this.f = this.d + Math.sqrt(this.d * this.d - 1.0);
4365          } else {
4366             this.f = this.d - Math.sqrt(this.d * this.d - 1.0);
4367          }
4368       } else {
4369         this.f = this.d;
4370       }
4371       this.el = this.f * Math.pow(this.ts,this.bl);
4372    }
4373
4374    //this.longc=52.60353916666667;
4375
4376    if (this.mode != 0) {
4377       this.g = .5 * (this.f - 1.0/this.f);
4378       this.gama = Proj4js.common.asinz(Math.sin(this.alpha) / this.d);
4379       this.longc= this.longc - Proj4js.common.asinz(this.g * Math.tan(this.gama))/this.bl;
4380
4381       /* Report parameters common to format B
4382       -------------------------------------*/
4383       //genrpt(azimuth * R2D,"Azimuth of Central Line:    ");
4384       //cenlon(lon_origin);
4385      // cenlat(lat_origin);
4386
4387       this.con = Math.abs(this.lat0);
4388       if ((this.con > Proj4js.common.EPSLN) && (Math.abs(this.con - Proj4js.common.HALF_PI) > Proj4js.common.EPSLN)) {
4389            this.singam=Math.sin(this.gama);
4390            this.cosgam=Math.cos(this.gama);
4391
4392            this.sinaz=Math.sin(this.alpha);
4393            this.cosaz=Math.cos(this.alpha);
4394
4395            if (this.lat0>= 0) {
4396               this.u =  (this.al / this.bl) * Math.atan(Math.sqrt(this.d*this.d - 1.0)/this.cosaz);
4397            } else {
4398               this.u =  -(this.al / this.bl) *Math.atan(Math.sqrt(this.d*this.d - 1.0)/this.cosaz);
4399            }
4400          } else {
4401            Proj4js.reportError("omerc:Init:DataError");
4402          }
4403       } else {
4404       this.sinphi =Math. sin(this.at1);
4405       this.ts1 = Proj4js.common.tsfnz(this.e,this.lat1,this.sinphi);
4406       this.sinphi = Math.sin(this.lat2);
4407       this.ts2 = Proj4js.common.tsfnz(this.e,this.lat2,this.sinphi);
4408       this.h = Math.pow(this.ts1,this.bl);
4409       this.l = Math.pow(this.ts2,this.bl);
4410       this.f = this.el/this.h;
4411       this.g = .5 * (this.f - 1.0/this.f);
4412       this.j = (this.el * this.el - this.l * this.h)/(this.el * this.el + this.l * this.h);
4413       this.p = (this.l - this.h) / (this.l + this.h);
4414       this.dlon = this.lon1 - this.lon2;
4415       if (this.dlon < -Proj4js.common.PI) this.lon2 = this.lon2 - 2.0 * Proj4js.common.PI;
4416       if (this.dlon > Proj4js.common.PI) this.lon2 = this.lon2 + 2.0 * Proj4js.common.PI;
4417       this.dlon = this.lon1 - this.lon2;
4418       this.longc = .5 * (this.lon1 + this.lon2) -Math.atan(this.j * Math.tan(.5 * this.bl * this.dlon)/this.p)/this.bl;
4419       this.dlon  = Proj4js.common.adjust_lon(this.lon1 - this.longc);
4420       this.gama = Math.atan(Math.sin(this.bl * this.dlon)/this.g);
4421       this.alpha = Proj4js.common.asinz(this.d * Math.sin(this.gama));
4422
4423       /* Report parameters common to format A
4424       -------------------------------------*/
4425
4426       if (Math.abs(this.lat1 - this.lat2) <= Proj4js.common.EPSLN) {
4427          Proj4js.reportError("omercInitDataError");
4428          //return(202);
4429       } else {
4430          this.con = Math.abs(this.lat1);
4431       }
4432       if ((this.con <= Proj4js.common.EPSLN) || (Math.abs(this.con - HALF_PI) <= Proj4js.common.EPSLN)) {
4433           Proj4js.reportError("omercInitDataError");
4434                //return(202);
4435       } else {
4436         if (Math.abs(Math.abs(this.lat0) - Proj4js.common.HALF_PI) <= Proj4js.common.EPSLN) {
4437            Proj4js.reportError("omercInitDataError");
4438            //return(202);
4439         }
4440       }
4441
4442       this.singam=Math.sin(this.gam);
4443       this.cosgam=Math.cos(this.gam);
4444
4445       this.sinaz=Math.sin(this.alpha);
4446       this.cosaz=Math.cos(this.alpha); 
4447
4448
4449       if (this.lat0 >= 0) {
4450          this.u =  (this.al/this.bl) * Math.atan(Math.sqrt(this.d * this.d - 1.0)/this.cosaz);
4451       } else {
4452          this.u = -(this.al/this.bl) * Math.atan(Math.sqrt(this.d * this.d - 1.0)/this.cosaz);
4453       }
4454     }
4455  },
4456
4457
4458  /* Oblique Mercator forward equations--mapping lat,long to x,y
4459    ----------------------------------------------------------*/
4460  forward: function(p) {
4461    var theta;          /* angle                                        */
4462    var sin_phi, cos_phi;/* sin and cos value                           */
4463    var b;              /* temporary values                             */
4464    var c, t, tq;       /* temporary values                             */
4465    var con, n, ml;     /* cone constant, small m                       */
4466    var q,us,vl;
4467    var ul,vs;
4468    var s;
4469    var dlon;
4470    var ts1;
4471
4472    var lon=p.x;
4473    var lat=p.y;
4474    /* Forward equations
4475      -----------------*/
4476    sin_phi = Math.sin(lat);
4477    dlon = Proj4js.common.adjust_lon(lon - this.longc);
4478    vl = Math.sin(this.bl * dlon);
4479    if (Math.abs(Math.abs(lat) - Proj4js.common.HALF_PI) > Proj4js.common.EPSLN) {
4480       ts1 = Proj4js.common.tsfnz(this.e,lat,sin_phi);
4481       q = this.el / (Math.pow(ts1,this.bl));
4482       s = .5 * (q - 1.0 / q);
4483       t = .5 * (q + 1.0/ q);
4484       ul = (s * this.singam - vl * this.cosgam) / t;
4485       con = Math.cos(this.bl * dlon);
4486       if (Math.abs(con) < .0000001) {
4487          us = this.al * this.bl * dlon;
4488       } else {
4489          us = this.al * Math.atan((s * this.cosgam + vl * this.singam) / con)/this.bl;
4490          if (con < 0) us = us + Proj4js.common.PI * this.al / this.bl;
4491       }
4492    } else {
4493       if (lat >= 0) {
4494          ul = this.singam;
4495       } else {
4496          ul = -this.singam;
4497       }
4498       us = this.al * lat / this.bl;
4499    }
4500    if (Math.abs(Math.abs(ul) - 1.0) <= Proj4js.common.EPSLN) {
4501       //alert("Point projects into infinity","omer-for");
4502       Proj4js.reportError("omercFwdInfinity");
4503       //return(205);
4504    }
4505    vs = .5 * this.al * Math.log((1.0 - ul)/(1.0 + ul)) / this.bl;
4506    us = us - this.u;
4507    var x = this.x0 + vs * this.cosaz + us * this.sinaz;
4508    var y = this.y0 + us * this.cosaz - vs * this.sinaz;
4509
4510    p.x=x;
4511    p.y=y;
4512    return p;
4513  },
4514
4515  inverse: function(p) {
4516    var delta_lon;      /* Delta longitude (Given longitude - center    */
4517    var theta;          /* angle                                        */
4518    var delta_theta;    /* adjusted longitude                           */
4519    var sin_phi, cos_phi;/* sin and cos value                           */
4520    var b;              /* temporary values                             */
4521    var c, t, tq;       /* temporary values                             */
4522    var con, n, ml;     /* cone constant, small m                       */
4523    var vs,us,q,s,ts1;
4524    var vl,ul,bs;
4525    var dlon;
4526    var  flag;
4527
4528    /* Inverse equations
4529      -----------------*/
4530    p.x -= this.x0;
4531    p.y -= this.y0;
4532    flag = 0;
4533    vs = p.x * this.cosaz - p.y * this.sinaz;
4534    us = p.y * this.cosaz + p.x * this.sinaz;
4535    us = us + this.u;
4536    q = Math.exp(-this.bl * vs / this.al);
4537    s = .5 * (q - 1.0/q);
4538    t = .5 * (q + 1.0/q);
4539    vl = Math.sin(this.bl * us / this.al);
4540    ul = (vl * this.cosgam + s * this.singam)/t;
4541    if (Math.abs(Math.abs(ul) - 1.0) <= Proj4js.common.EPSLN)
4542       {
4543       lon = this.longc;
4544       if (ul >= 0.0) {
4545          lat = Proj4js.common.HALF_PI;
4546       } else {
4547         lat = -Proj4js.common.HALF_PI;
4548       }
4549    } else {
4550       con = 1.0 / this.bl;
4551       ts1 =Math.pow((this.el / Math.sqrt((1.0 + ul) / (1.0 - ul))),con);
4552       lat = Proj4js.common.phi2z(this.e,ts1);
4553       //if (flag != 0)
4554          //return(flag);
4555       //~ con = Math.cos(this.bl * us /al);
4556       theta = this.longc - Math.atan2((s * this.cosgam - vl * this.singam) , con)/this.bl;
4557       lon = Proj4js.common.adjust_lon(theta);
4558    }
4559    p.x=lon;
4560    p.y=lat;
4561    return p;
4562  }
4563};
4564/* ======================================================================
4565    projCode/lcc.js
4566   ====================================================================== */
4567
4568/*******************************************************************************
4569NAME                            LAMBERT CONFORMAL CONIC
4570
4571PURPOSE:        Transforms input longitude and latitude to Easting and
4572                Northing for the Lambert Conformal Conic projection.  The
4573                longitude and latitude must be in radians.  The Easting
4574                and Northing values will be returned in meters.
4575
4576
4577ALGORITHM REFERENCES
4578
45791.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
4580    Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
4581    State Government Printing Office, Washington D.C., 1987.
4582
45832.  Snyder, John P. and Voxland, Philip M., "An Album of Map Projections",
4584    U.S. Geological Survey Professional Paper 1453 , United State Government
4585*******************************************************************************/
4586
4587
4588//<2104> +proj=lcc +lat_1=10.16666666666667 +lat_0=10.16666666666667 +lon_0=-71.60561777777777 +k_0=1 +x0=-17044 +x0=-23139.97 +ellps=intl +units=m +no_defs  no_defs
4589
4590// Initialize the Lambert Conformal conic projection
4591// -----------------------------------------------------------------
4592
4593//Proj4js.Proj.lcc = Class.create();
4594Proj4js.Proj.lcc = {
4595  init : function() {
4596
4597    // array of:  r_maj,r_min,lat1,lat2,c_lon,c_lat,false_east,false_north
4598    //double c_lat;                   /* center latitude                      */
4599    //double c_lon;                   /* center longitude                     */
4600    //double lat1;                    /* first standard parallel              */
4601    //double lat2;                    /* second standard parallel             */
4602    //double r_maj;                   /* major axis                           */
4603    //double r_min;                   /* minor axis                           */
4604    //double false_east;              /* x offset in meters                   */
4605    //double false_north;             /* y offset in meters                   */
4606
4607      if (!this.lat2){this.lat2=this.lat0;}//if lat2 is not defined
4608      if (!this.k0) this.k0 = 1.0;
4609
4610    // Standard Parallels cannot be equal and on opposite sides of the equator
4611      if (Math.abs(this.lat1+this.lat2) < Proj4js.common.EPSLN) {
4612        Proj4js.reportError("lcc:init: Equal Latitudes");
4613        return;
4614      }
4615
4616      var temp = this.b / this.a;
4617      this.e = Math.sqrt(1.0 - temp*temp);
4618
4619      var sin1 = Math.sin(this.lat1);
4620      var cos1 = Math.cos(this.lat1);
4621      var ms1 = Proj4js.common.msfnz(this.e, sin1, cos1);
4622      var ts1 = Proj4js.common.tsfnz(this.e, this.lat1, sin1);
4623
4624      var sin2 = Math.sin(this.lat2);
4625      var cos2 = Math.cos(this.lat2);
4626      var ms2 = Proj4js.common.msfnz(this.e, sin2, cos2);
4627      var ts2 = Proj4js.common.tsfnz(this.e, this.lat2, sin2);
4628
4629      var ts0 = Proj4js.common.tsfnz(this.e, this.lat0, Math.sin(this.lat0));
4630
4631      if (Math.abs(this.lat1 - this.lat2) > Proj4js.common.EPSLN) {
4632        this.ns = Math.log(ms1/ms2)/Math.log(ts1/ts2);
4633      } else {
4634        this.ns = sin1;
4635      }
4636      this.f0 = ms1 / (this.ns * Math.pow(ts1, this.ns));
4637      this.rh = this.a * this.f0 * Math.pow(ts0, this.ns);
4638      if (!this.title) this.title = "Lambert Conformal Conic";
4639    },
4640
4641
4642    // Lambert Conformal conic forward equations--mapping lat,long to x,y
4643    // -----------------------------------------------------------------
4644    forward : function(p) {
4645
4646      var lon = p.x;
4647      var lat = p.y;
4648
4649    // convert to radians
4650      if ( lat <= 90.0 && lat >= -90.0 && lon <= 180.0 && lon >= -180.0) {
4651        //lon = lon * Proj4js.common.D2R;
4652        //lat = lat * Proj4js.common.D2R;
4653      } else {
4654        Proj4js.reportError("lcc:forward: llInputOutOfRange: "+ lon +" : " + lat);
4655        return null;
4656      }
4657
4658      var con  = Math.abs( Math.abs(lat) - Proj4js.common.HALF_PI);
4659      var ts, rh1;
4660      if (con > Proj4js.common.EPSLN) {
4661        ts = Proj4js.common.tsfnz(this.e, lat, Math.sin(lat) );
4662        rh1 = this.a * this.f0 * Math.pow(ts, this.ns);
4663      } else {
4664        con = lat * this.ns;
4665        if (con <= 0) {
4666          Proj4js.reportError("lcc:forward: No Projection");
4667          return null;
4668        }
4669        rh1 = 0;
4670      }
4671      var theta = this.ns * Proj4js.common.adjust_lon(lon - this.long0);
4672      p.x = this.k0 * (rh1 * Math.sin(theta)) + this.x0;
4673      p.y = this.k0 * (this.rh - rh1 * Math.cos(theta)) + this.y0;
4674
4675      return p;
4676    },
4677
4678  // Lambert Conformal Conic inverse equations--mapping x,y to lat/long
4679  // -----------------------------------------------------------------
4680  inverse : function(p) {
4681
4682    var rh1, con, ts;
4683    var lat, lon;
4684    var x = (p.x - this.x0)/this.k0;
4685    var y = (this.rh - (p.y - this.y0)/this.k0);
4686    if (this.ns > 0) {
4687      rh1 = Math.sqrt (x * x + y * y);
4688      con = 1.0;
4689    } else {
4690      rh1 = -Math.sqrt (x * x + y * y);
4691      con = -1.0;
4692    }
4693    var theta = 0.0;
4694    if (rh1 != 0) {
4695      theta = Math.atan2((con * x),(con * y));
4696    }
4697    if ((rh1 != 0) || (this.ns > 0.0)) {
4698      con = 1.0/this.ns;
4699      ts = Math.pow((rh1/(this.a * this.f0)), con);
4700      lat = Proj4js.common.phi2z(this.e, ts);
4701      if (lat == -9999) return null;
4702    } else {
4703      lat = -Proj4js.common.HALF_PI;
4704    }
4705    lon = Proj4js.common.adjust_lon(theta/this.ns + this.long0);
4706
4707    p.x = lon;
4708    p.y = lat;
4709    return p;
4710  }
4711};
4712
4713
4714
4715
4716/* ======================================================================
4717    projCode/laea.js
4718   ====================================================================== */
4719
4720/*******************************************************************************
4721NAME                  LAMBERT AZIMUTHAL EQUAL-AREA
4722 
4723PURPOSE:        Transforms input longitude and latitude to Easting and
4724                Northing for the Lambert Azimuthal Equal-Area projection.  The
4725                longitude and latitude must be in radians.  The Easting
4726                and Northing values will be returned in meters.
4727
4728PROGRAMMER              DATE           
4729----------              ----           
4730D. Steinwand, EROS      March, 1991   
4731
4732This function was adapted from the Lambert Azimuthal Equal Area projection
4733code (FORTRAN) in the General Cartographic Transformation Package software
4734which is available from the U.S. Geological Survey National Mapping Division.
4735 
4736ALGORITHM REFERENCES
4737
47381.  "New Equal-Area Map Projections for Noncircular Regions", John P. Snyder,
4739    The American Cartographer, Vol 15, No. 4, October 1988, pp. 341-355.
4740
47412.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
4742    Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
4743    State Government Printing Office, Washington D.C., 1987.
4744
47453.  "Software Documentation for GCTP General Cartographic Transformation
4746    Package", U.S. Geological Survey National Mapping Division, May 1982.
4747*******************************************************************************/
4748
4749Proj4js.Proj.laea = {
4750  S_POLE: 1,
4751  N_POLE: 2,
4752  EQUIT: 3,
4753  OBLIQ: 4,
4754
4755
4756/* Initialize the Lambert Azimuthal Equal Area projection
4757  ------------------------------------------------------*/
4758  init: function() {
4759    var t = Math.abs(this.lat0);
4760    if (Math.abs(t - Proj4js.common.HALF_PI) < Proj4js.common.EPSLN) {
4761      this.mode = this.lat0 < 0. ? this.S_POLE : this.N_POLE;
4762    } else if (Math.abs(t) < Proj4js.common.EPSLN) {
4763      this.mode = this.EQUIT;
4764    } else {
4765      this.mode = this.OBLIQ;
4766    }
4767    if (this.es > 0) {
4768      var sinphi;
4769 
4770      this.qp = Proj4js.common.qsfnz(this.e, 1.0);
4771      this.mmf = .5 / (1. - this.es);
4772      this.apa = this.authset(this.es);
4773      switch (this.mode) {
4774        case this.N_POLE:
4775        case this.S_POLE:
4776          this.dd = 1.;
4777          break;
4778        case this.EQUIT:
4779          this.rq = Math.sqrt(.5 * this.qp);
4780          this.dd = 1. / this.rq;
4781          this.xmf = 1.;
4782          this.ymf = .5 * this.qp;
4783          break;
4784        case this.OBLIQ:
4785          this.rq = Math.sqrt(.5 * this.qp);
4786          sinphi = Math.sin(this.lat0);
4787          this.sinb1 = Proj4js.common.qsfnz(this.e, sinphi) / this.qp;
4788          this.cosb1 = Math.sqrt(1. - this.sinb1 * this.sinb1);
4789          this.dd = Math.cos(this.lat0) / (Math.sqrt(1. - this.es * sinphi * sinphi) * this.rq * this.cosb1);
4790          this.ymf = (this.xmf = this.rq) / this.dd;
4791          this.xmf *= this.dd;
4792          break;
4793      }
4794    } else {
4795      if (this.mode == this.OBLIQ) {
4796        this.sinph0 = Math.sin(this.lat0);
4797        this.cosph0 = Math.cos(this.lat0);
4798      }
4799    }
4800  },
4801
4802/* Lambert Azimuthal Equal Area forward equations--mapping lat,long to x,y
4803  -----------------------------------------------------------------------*/
4804  forward: function(p) {
4805
4806    /* Forward equations
4807      -----------------*/
4808    var x,y;
4809    var lam=p.x;
4810    var phi=p.y;
4811    lam = Proj4js.common.adjust_lon(lam - this.long0);
4812   
4813    if (this.sphere) {
4814        var coslam, cosphi, sinphi;
4815     
4816        sinphi = Math.sin(phi);
4817        cosphi = Math.cos(phi);
4818        coslam = Math.cos(lam);
4819        switch (this.mode) {
4820          case this.OBLIQ:
4821          case this.EQUIT:
4822            y = (this.mode == this.EQUIT) ? 1. + cosphi * coslam : 1. + this.sinph0 * sinphi + this.cosph0 * cosphi * coslam;
4823            if (y <= Proj4js.common.EPSLN) {
4824              Proj4js.reportError("laea:fwd:y less than eps");
4825              return null;
4826            }
4827            y = Math.sqrt(2. / y);
4828            x = y * cosphi * Math.sin(lam);
4829            y *= (this.mode == this.EQUIT) ? sinphi : this.cosph0 * sinphi - this.sinph0 * cosphi * coslam;
4830            break;
4831          case this.N_POLE:
4832            coslam = -coslam;
4833          case this.S_POLE:
4834            if (Math.abs(phi + this.phi0) < Proj4js.common.EPSLN) {
4835              Proj4js.reportError("laea:fwd:phi < eps");
4836              return null;
4837            }
4838            y = Proj4js.common.FORTPI - phi * .5;
4839            y = 2. * ((this.mode == this.S_POLE) ? Math.cos(y) : Math.sin(y));
4840            x = y * Math.sin(lam);
4841            y *= coslam;
4842            break;
4843        }
4844    } else {
4845        var coslam, sinlam, sinphi, q, sinb=0.0, cosb=0.0, b=0.0;
4846     
4847        coslam = Math.cos(lam);
4848        sinlam = Math.sin(lam);
4849        sinphi = Math.sin(phi);
4850        q = Proj4js.common.qsfnz(this.e, sinphi);
4851        if (this.mode == this.OBLIQ || this.mode == this.EQUIT) {
4852          sinb = q / this.qp;
4853          cosb = Math.sqrt(1. - sinb * sinb);
4854        }
4855        switch (this.mode) {
4856          case this.OBLIQ:
4857            b = 1. + this.sinb1 * sinb + this.cosb1 * cosb * coslam;
4858            break;
4859          case this.EQUIT:
4860            b = 1. + cosb * coslam;
4861            break;
4862          case this.N_POLE:
4863            b = Proj4js.common.HALF_PI + phi;
4864            q = this.qp - q;
4865            break;
4866          case this.S_POLE:
4867            b = phi - Proj4js.common.HALF_PI;
4868            q = this.qp + q;
4869            break;
4870        }
4871        if (Math.abs(b) < Proj4js.common.EPSLN) {
4872            Proj4js.reportError("laea:fwd:b < eps");
4873            return null;
4874        }
4875        switch (this.mode) {
4876          case this.OBLIQ:
4877          case this.EQUIT:
4878            b = Math.sqrt(2. / b);
4879            if (this.mode == this.OBLIQ) {
4880              y = this.ymf * b * (this.cosb1 * sinb - this.sinb1 * cosb * coslam);
4881            } else {
4882              y = (b = Math.sqrt(2. / (1. + cosb * coslam))) * sinb * this.ymf;
4883            }
4884            x = this.xmf * b * cosb * sinlam;
4885            break;
4886          case this.N_POLE:
4887          case this.S_POLE:
4888            if (q >= 0.) {
4889              x = (b = Math.sqrt(q)) * sinlam;
4890              y = coslam * ((this.mode == this.S_POLE) ? b : -b);
4891            } else {
4892              x = y = 0.;
4893            }
4894            break;
4895        }
4896    }
4897
4898    //v 1.0
4899    /*
4900    var sin_lat=Math.sin(lat);
4901    var cos_lat=Math.cos(lat);
4902
4903    var sin_delta_lon=Math.sin(delta_lon);
4904    var cos_delta_lon=Math.cos(delta_lon);
4905
4906    var g =this.sin_lat_o * sin_lat +this.cos_lat_o * cos_lat * cos_delta_lon;
4907    if (g == -1.0) {
4908      Proj4js.reportError("laea:fwd:Point projects to a circle of radius "+ 2.0 * R);
4909      return null;
4910    }
4911    var ksp = this.a * Math.sqrt(2.0 / (1.0 + g));
4912    var x = ksp * cos_lat * sin_delta_lon + this.x0;
4913    var y = ksp * (this.cos_lat_o * sin_lat - this.sin_lat_o * cos_lat * cos_delta_lon) + this.y0;
4914    */
4915    p.x = this.a*x + this.x0;
4916    p.y = this.a*y + this.y0;
4917    return p;
4918  },//lamazFwd()
4919
4920/* Inverse equations
4921  -----------------*/
4922  inverse: function(p) {
4923    p.x -= this.x0;
4924    p.y -= this.y0;
4925    var x = p.x/this.a;
4926    var y = p.y/this.a;
4927   
4928    if (this.sphere) {
4929        var  cosz=0.0, rh, sinz=0.0;
4930     
4931        rh = Math.sqrt(x*x + y*y);
4932        var phi = rh * .5;
4933        if (phi > 1.) {
4934          Proj4js.reportError("laea:Inv:DataError");
4935          return null;
4936        }
4937        phi = 2. * Math.asin(phi);
4938        if (this.mode == this.OBLIQ || this.mode == this.EQUIT) {
4939          sinz = Math.sin(phi);
4940          cosz = Math.cos(phi);
4941        }
4942        switch (this.mode) {
4943        case this.EQUIT:
4944          phi = (Math.abs(rh) <= Proj4js.common.EPSLN) ? 0. : Math.asin(y * sinz / rh);
4945          x *= sinz;
4946          y = cosz * rh;
4947          break;
4948        case this.OBLIQ:
4949          phi = (Math.abs(rh) <= Proj4js.common.EPSLN) ? this.phi0 : Math.asin(cosz * sinph0 + y * sinz * cosph0 / rh);
4950          x *= sinz * cosph0;
4951          y = (cosz - Math.sin(phi) * sinph0) * rh;
4952          break;
4953        case this.N_POLE:
4954          y = -y;
4955          phi = Proj4js.common.HALF_PI - phi;
4956          break;
4957        case this.S_POLE:
4958          phi -= Proj4js.common.HALF_PI;
4959          break;
4960        }
4961        lam = (y == 0. && (this.mode == this.EQUIT || this.mode == this.OBLIQ)) ? 0. : Math.atan2(x, y);
4962    } else {
4963        var cCe, sCe, q, rho, ab=0.0;
4964     
4965        switch (this.mode) {
4966          case this.EQUIT:
4967          case this.OBLIQ:
4968            x /= this.dd;
4969            y *=  this.dd;
4970            rho = Math.sqrt(x*x + y*y);
4971            if (rho < Proj4js.common.EPSLN) {
4972              p.x = 0.;
4973              p.y = this.phi0;
4974              return p;
4975            }
4976            sCe = 2. * Math.asin(.5 * rho / this.rq);
4977            cCe = Math.cos(sCe);
4978            x *= (sCe = Math.sin(sCe));
4979            if (this.mode == this.OBLIQ) {
4980              ab = cCe * this.sinb1 + y * sCe * this.cosb1 / rho
4981              q = this.qp * ab;
4982              y = rho * this.cosb1 * cCe - y * this.sinb1 * sCe;
4983            } else {
4984              ab = y * sCe / rho;
4985              q = this.qp * ab;
4986              y = rho * cCe;
4987            }
4988            break;
4989          case this.N_POLE:
4990            y = -y;
4991          case this.S_POLE:
4992            q = (x * x + y * y);
4993            if (!q ) {
4994              p.x = 0.;
4995              p.y = this.phi0;
4996              return p;
4997            }
4998            /*
4999            q = this.qp - q;
5000            */
5001            ab = 1. - q / this.qp;
5002            if (this.mode == this.S_POLE) {
5003              ab = - ab;
5004            }
5005            break;
5006        }
5007        lam = Math.atan2(x, y);
5008        phi = this.authlat(Math.asin(ab), this.apa);
5009    }
5010
5011    /*
5012    var Rh = Math.Math.sqrt(p.x *p.x +p.y * p.y);
5013    var temp = Rh / (2.0 * this.a);
5014
5015    if (temp > 1) {
5016      Proj4js.reportError("laea:Inv:DataError");
5017      return null;
5018    }
5019
5020    var z = 2.0 * Proj4js.common.asinz(temp);
5021    var sin_z=Math.sin(z);
5022    var cos_z=Math.cos(z);
5023
5024    var lon =this.long0;
5025    if (Math.abs(Rh) > Proj4js.common.EPSLN) {
5026       var lat = Proj4js.common.asinz(this.sin_lat_o * cos_z +this. cos_lat_o * sin_z *p.y / Rh);
5027       var temp =Math.abs(this.lat0) - Proj4js.common.HALF_PI;
5028       if (Math.abs(temp) > Proj4js.common.EPSLN) {
5029          temp = cos_z -this.sin_lat_o * Math.sin(lat);
5030          if(temp!=0.0) lon=Proj4js.common.adjust_lon(this.long0+Math.atan2(p.x*sin_z*this.cos_lat_o,temp*Rh));
5031       } else if (this.lat0 < 0.0) {
5032          lon = Proj4js.common.adjust_lon(this.long0 - Math.atan2(-p.x,p.y));
5033       } else {
5034          lon = Proj4js.common.adjust_lon(this.long0 + Math.atan2(p.x, -p.y));
5035       }
5036    } else {
5037      lat = this.lat0;
5038    }
5039    */
5040    //return(OK);
5041    p.x = Proj4js.common.adjust_lon(this.long0+lam);
5042    p.y = phi;
5043    return p;
5044  },//lamazInv()
5045 
5046/* determine latitude from authalic latitude */
5047  P00: .33333333333333333333,
5048  P01: .17222222222222222222,
5049  P02: .10257936507936507936,
5050  P10: .06388888888888888888,
5051  P11: .06640211640211640211,
5052  P20: .01641501294219154443,
5053 
5054  authset: function(es) {
5055    var t;
5056    var APA = new Array();
5057    APA[0] = es * this.P00;
5058    t = es * es;
5059    APA[0] += t * this.P01;
5060    APA[1] = t * this.P10;
5061    t *= es;
5062    APA[0] += t * this.P02;
5063    APA[1] += t * this.P11;
5064    APA[2] = t * this.P20;
5065    return APA;
5066  },
5067 
5068  authlat: function(beta, APA) {
5069    var t = beta+beta;
5070    return(beta + APA[0] * Math.sin(t) + APA[1] * Math.sin(t+t) + APA[2] * Math.sin(t+t+t));
5071  }
5072 
5073};
5074
5075
5076
5077/* ======================================================================
5078    projCode/aeqd.js
5079   ====================================================================== */
5080
5081Proj4js.Proj.aeqd = {
5082
5083  init : function() {
5084    this.sin_p12=Math.sin(this.lat0);
5085    this.cos_p12=Math.cos(this.lat0);
5086  },
5087
5088  forward: function(p) {
5089    var lon=p.x;
5090    var lat=p.y;
5091    var ksp;
5092
5093    var sinphi=Math.sin(p.y);
5094    var cosphi=Math.cos(p.y); 
5095    var dlon = Proj4js.common.adjust_lon(lon - this.long0);
5096    var coslon = Math.cos(dlon);
5097    var g = this.sin_p12 * sinphi + this.cos_p12 * cosphi * coslon;
5098    if (Math.abs(Math.abs(g) - 1.0) < Proj4js.common.EPSLN) {
5099       ksp = 1.0;
5100       if (g < 0.0) {
5101         Proj4js.reportError("aeqd:Fwd:PointError");
5102         return;
5103       }
5104    } else {
5105       var z = Math.acos(g);
5106       ksp = z/Math.sin(z);
5107    }
5108    p.x = this.x0 + this.a * ksp * cosphi * Math.sin(dlon);
5109    p.y = this.y0 + this.a * ksp * (this.cos_p12 * sinphi - this.sin_p12 * cosphi * coslon);
5110    return p;
5111  },
5112
5113  inverse: function(p){
5114    p.x -= this.x0;
5115    p.y -= this.y0;
5116
5117    var rh = Math.sqrt(p.x * p.x + p.y *p.y);
5118    if (rh > (2.0 * Proj4js.common.HALF_PI * this.a)) {
5119       Proj4js.reportError("aeqdInvDataError");
5120       return;
5121    }
5122    var z = rh / this.a;
5123
5124    var sinz=Math.sin(z);
5125    var cosz=Math.cos(z);
5126
5127    var lon = this.long0;
5128    var lat;
5129    if (Math.abs(rh) <= Proj4js.common.EPSLN) {
5130      lat = this.lat0;
5131    } else {
5132      lat = Proj4js.common.asinz(cosz * this.sin_p12 + (p.y * sinz * this.cos_p12) / rh);
5133      var con = Math.abs(this.lat0) - Proj4js.common.HALF_PI;
5134      if (Math.abs(con) <= Proj4js.common.EPSLN) {
5135        if (lat0 >= 0.0) {
5136          lon = Proj4js.common.adjust_lon(this.long0 + Math.atan2(p.x , -p.y));
5137        } else {
5138          lon = Proj4js.common.adjust_lon(this.long0 - Math.atan2(-p.x , p.y));
5139        }
5140      } else {
5141        con = cosz - this.sin_p12 * Math.sin(lat);
5142        if ((Math.abs(con) < Proj4js.common.EPSLN) && (Math.abs(p.x) < Proj4js.common.EPSLN)) {
5143           //no-op, just keep the lon value as is
5144        } else {
5145          var temp = Math.atan2((p.x * sinz * this.cos_p12), (con * rh));
5146          lon = Proj4js.common.adjust_lon(this.long0 + Math.atan2((p.x * sinz * this.cos_p12), (con * rh)));
5147        }
5148      }
5149    }
5150
5151    p.x = lon;
5152    p.y = lat;
5153    return p;
5154  } 
5155};
5156/* ======================================================================
5157    projCode/moll.js
5158   ====================================================================== */
5159
5160/*******************************************************************************
5161NAME                            MOLLWEIDE
5162
5163PURPOSE:        Transforms input longitude and latitude to Easting and
5164                Northing for the MOllweide projection.  The
5165                longitude and latitude must be in radians.  The Easting
5166                and Northing values will be returned in meters.
5167
5168PROGRAMMER              DATE
5169----------              ----
5170D. Steinwand, EROS      May, 1991;  Updated Sept, 1992; Updated Feb, 1993
5171S. Nelson, EDC          Jun, 2993;      Made corrections in precision and
5172                                        number of iterations.
5173
5174ALGORITHM REFERENCES
5175
51761.  Snyder, John P. and Voxland, Philip M., "An Album of Map Projections",
5177    U.S. Geological Survey Professional Paper 1453 , United State Government
5178    Printing Office, Washington D.C., 1989.
5179
51802.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
5181    Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
5182    State Government Printing Office, Washington D.C., 1987.
5183*******************************************************************************/
5184
5185Proj4js.Proj.moll = {
5186
5187  /* Initialize the Mollweide projection
5188    ------------------------------------*/
5189  init: function(){
5190    //no-op
5191  },
5192
5193  /* Mollweide forward equations--mapping lat,long to x,y
5194    ----------------------------------------------------*/
5195  forward: function(p) {
5196
5197    /* Forward equations
5198      -----------------*/
5199    var lon=p.x;
5200    var lat=p.y;
5201
5202    var delta_lon = Proj4js.common.adjust_lon(lon - this.long0);
5203    var theta = lat;
5204    var con = Proj4js.common.PI * Math.sin(lat);
5205
5206    /* Iterate using the Newton-Raphson method to find theta
5207      -----------------------------------------------------*/
5208    for (var i=0;true;i++) {
5209       var delta_theta = -(theta + Math.sin(theta) - con)/ (1.0 + Math.cos(theta));
5210       theta += delta_theta;
5211       if (Math.abs(delta_theta) < Proj4js.common.EPSLN) break;
5212       if (i >= 50) {
5213          Proj4js.reportError("moll:Fwd:IterationError");
5214         //return(241);
5215       }
5216    }
5217    theta /= 2.0;
5218
5219    /* If the latitude is 90 deg, force the x coordinate to be "0 + false easting"
5220       this is done here because of precision problems with "cos(theta)"
5221       --------------------------------------------------------------------------*/
5222    if (Proj4js.common.PI/2 - Math.abs(lat) < Proj4js.common.EPSLN) delta_lon =0;
5223    var x = 0.900316316158 * this.a * delta_lon * Math.cos(theta) + this.x0;
5224    var y = 1.4142135623731 * this.a * Math.sin(theta) + this.y0;
5225
5226    p.x=x;
5227    p.y=y;
5228    return p;
5229  },
5230
5231  inverse: function(p){
5232    var theta;
5233    var arg;
5234
5235    /* Inverse equations
5236      -----------------*/
5237    p.x-= this.x0;
5238    //~ p.y -= this.y0;
5239    var arg = p.y /  (1.4142135623731 * this.a);
5240
5241    /* Because of division by zero problems, 'arg' can not be 1.0.  Therefore
5242       a number very close to one is used instead.
5243       -------------------------------------------------------------------*/
5244    if(Math.abs(arg) > 0.999999999999) arg=0.999999999999;
5245    var theta =Math.asin(arg);
5246    var lon = Proj4js.common.adjust_lon(this.long0 + (p.x / (0.900316316158 * this.a * Math.cos(theta))));
5247    if(lon < (-Proj4js.common.PI)) lon= -Proj4js.common.PI;
5248    if(lon > Proj4js.common.PI) lon= Proj4js.common.PI;
5249    arg = (2.0 * theta + Math.sin(2.0 * theta)) / Proj4js.common.PI;
5250    if(Math.abs(arg) > 1.0)arg=1.0;
5251    var lat = Math.asin(arg);
5252    //return(OK);
5253
5254    p.x=lon;
5255    p.y=lat;
5256    return p;
5257  }
5258};
5259
Note: See TracBrowser for help on using the repository browser.