001/*
002 * To change this template, choose Tools | Templates
003 * and open the template in the editor.
004 */
005
006package armyc2.c2sd.JavaTacticalRenderer;
007import armyc2.c2sd.JavaLineArray.POINT2;
008import armyc2.c2sd.JavaLineArray.ref;
009import java.util.ArrayList;
010import armyc2.c2sd.renderer.utilities.ErrorLogger;
011import armyc2.c2sd.renderer.utilities.RendererException;
012import armyc2.c2sd.graphics2d.Rectangle2D;
013/**
014 * Class to calculate the geodesic based shapes for the Fire Support Areas
015 * @author Michael Deutch
016 */
017public final class mdlGeodesic {
018    private static final String _className = "mdlGeodesic";
019    private static final double sm_a    = 6378137;
020
021    private static double DegToRad(double deg) {
022        return deg / 180.0 * Math.PI;
023    }
024
025    private static double RadToDeg(double rad) {
026        return rad / Math.PI * 180.0;
027    }
028/**
029 * Returns the azimuth from true north between two points
030 * @param c1
031 * @param c2
032 * @return the azimuth from c1 to c2
033 */
034    public static double GetAzimuth(POINT2 c1, 
035            POINT2 c2) {//was private
036        double theta = 0;
037        try {
038            double lat1 = DegToRad(c1.y);
039            double lon1 = DegToRad(c1.x);
040            double lat2 = DegToRad(c2.y);
041            double lon2 = DegToRad(c2.x);
042            //formula
043            //θ = atan2( sin(Δlong).cos(lat2),
044            //cos(lat1).sin(lat2) − sin(lat1).cos(lat2).cos(Δlong) )
045            //var theta:Number = Math.atan2( Math.sin(lon2-lon1)*Math.cos(lat2),
046            //Math.cos(lat1)*Math.sin(lat2) − Math.sin(lat1)*Math.cos(lat2)*Math.cos(lon2-lon1) );
047            double y = Math.sin(lon2 - lon1);
048            y *= Math.cos(lat2);
049            double x = Math.cos(lat1);
050            x *= Math.sin(lat2);
051            double z = Math.sin(lat1);
052            z *= Math.cos(lat2);
053            z *= Math.cos(lon2 - lon1);
054            x = x - z;
055            theta = Math.atan2(y, x);
056            theta = RadToDeg(theta);
057        }
058        catch (Exception exc) {
059            //System.out.println(e.getMessage());
060            //clsUtility.WriteFile("Error in mdlGeodesic.GetAzimuth");
061               ErrorLogger.LogException(_className ,"GetAzimuth",
062                    new RendererException("Failed inside GetAzimuth", exc));
063        }
064        return theta;//RadToDeg(k);
065    }
066    /**
067     * Calculates the distance in meters between two geodesic points.
068     * Also calculates the azimuth from c1 to c2 and from c2 to c1.
069     *
070     * @param c1 the first point
071     * @param c2 the last point
072     * @param a12 OUT - an object with a member to hold the calculated azimuth in degrees from c1 to c2
073     * @param a21 OUT - an object with a member to hold the calculated azimuth in degrees from c2 to c1
074     * @return the distance in meters between c1 and c2
075     */
076    public static double geodesic_distance(POINT2 c1,
077            POINT2 c2,
078            ref<double[]> a12,
079            ref<double[]> a21) {
080        double h = 0;
081        try {
082            //formula
083            //R = earth’s radius (mean radius = 6,371km)
084            //Δlat = lat2− lat1
085            //Δlong = long2− long1
086            //a = sin²(Δlat/2) + cos(lat1).cos(lat2).sin²(Δlong/2)
087            //c = 2.atan2(√a, √(1−a))
088            //d = R.c
089            if(a12 != null && a21 !=null)
090            {
091                a12.value = new double[1];
092                a21.value = new double[1];
093                //set the azimuth
094                a12.value[0] = GetAzimuth(c1, c2);
095                a21.value[0] = GetAzimuth(c2, c1);
096            }
097            //c1.x+=360;
098            double dLat = DegToRad(c2.y - c1.y);
099            double dLon = DegToRad(c2.x - c1.x);
100
101            double b = 0, lat1 = 0, lat2 = 0, e = 0, f = 0, g = 0, k = 0;
102            b = Math.sin(dLat / 2);
103            lat1 = DegToRad(c1.y);
104            lat2 = DegToRad(c2.y);
105            e = Math.sin(dLon / 2);
106            f = Math.cos(lat1);
107            g = Math.cos(lat2);
108            //uncomment this to test calculation
109            //var a:Number = Math.sin(dLat / 2) * Math.sin(dLat / 2) + Math.cos(DegToRad(c1.y)) * Math.cos(DegToRad(c2.y)) * Math.sin(dLon / 2) * Math.sin(dLon / 2);
110            double a = b * b + f * g * e * e;
111            h = Math.sqrt(a);
112            k = Math.sqrt(1 - a);
113            h = 2 * Math.atan2(h, k);
114        }
115        catch (Exception exc) {
116            //System.out.println(e.getMessage());
117            //clsUtility.WriteFile("Error in mdlGeodesic.geodesic_distance");
118               ErrorLogger.LogException(_className ,"geodesic_distance",
119                    new RendererException("Failed inside geodesic_distance", exc));
120        }
121        return sm_a * h;
122    }
123    /**
124     * Calculates a geodesic point and given distance and azimuth from the srating geodesic point
125     *
126     * @param start the starting point
127     * @param distance the distance in meters
128     * @param azimuth the azimuth or bearing in degrees
129     *
130     * @return the calculated point
131     */
132    public static POINT2 geodesic_coordinate(POINT2 start,
133            double distance,
134            double azimuth) {
135        POINT2 pt = null;
136        try
137        {
138        //formula
139        //lat2 = asin(sin(lat1)*cos(d/R) + cos(lat1)*sin(d/R)*cos(θ))
140        //lon2 = lon1 + atan2(sin(θ)*sin(d/R)*cos(lat1), cos(d/R)−sin(lat1)*sin(lat2))
141
142        double a = 0, b = 0, c = 0, d = 0, e = 0, f = 0, g = 0, h = 0,
143                j = 0, k = 0, l = 0, m = 0, n = 0, p = 0, q = 0;
144
145        a = DegToRad(start.y);
146        b = Math.cos(a);
147        c = DegToRad(azimuth);
148        d = Math.sin(a);
149        e = Math.cos(distance / sm_a);
150        f = Math.sin(distance / sm_a);
151        g = Math.cos(c);
152        //uncomment to test calculation
153        //var lat2:Number = RadToDeg(Math.asin(Math.sin(DegToRad(start.y)) * Math.cos(DegToRad(distance / sm_a)) + Math.cos(DegToRad(start.y)) * Math.sin(DegToRad(distance / sm_a)) * Math.cos(DegToRad(azimuth))));
154        //lat2 = asin(sin(lat1)*cos(d/R) + cos(lat1)*sin(d/R)*cos(θ))
155        //var lat2:Number = RadToDeg(Math.asin(Math.sin(DegToRad(start.y)) * Math.cos(distance / sm_a) + Math.cos(DegToRad(start.y)) * Math.sin(distance / sm_a) * Math.cos(DegToRad(azimuth))));
156        //double lat2 = RadToDeg(Math.asin(Math.sin(DegToRad(start.y)) * Math.cos(distance / sm_a) + Math.cos(DegToRad(start.y)) * Math.sin(distance / sm_a) * Math.cos(DegToRad(azimuth))));
157        double lat = RadToDeg(Math.asin(d * e + b * f * g));
158        h = Math.sin(c);
159        k = Math.sin(h);
160        l = Math.cos(a);
161        m = DegToRad(lat);
162        n = Math.sin(m);
163        p = Math.atan2(h * f * b, e - d * n);
164        //uncomment to test calculation
165        //var lon2:Number = start.x + DegToRad(Math.atan2(Math.sin(DegToRad(azimuth)) * Math.sin(DegToRad(distance / sm_a)) * Math.cos(DegToRad(start.y)), Math.cos(DegToRad(distance / sm_a)) - Math.sin(DegToRad(start.y)) * Math.sin(DegToRad(lat))));
166        //lon2 = lon1 + atan2(sin(θ)*sin(d/R)*cos(lat1), cos(d/R)−sin(lat1)*sin(lat2))
167        //var lon2:Number = start.x + RadToDeg(Math.atan2(Math.sin(DegToRad(azimuth)) * Math.sin(distance / sm_a) * Math.cos(DegToRad(start.y)), Math.cos(distance / sm_a) - Math.sin(DegToRad(start.y)) * Math.sin(DegToRad(lat2))));
168        double lon = start.x + RadToDeg(p);
169        pt = new POINT2(lon, lat);
170        }
171        catch (Exception exc) {
172            //clsUtility.WriteFile("Error in mdlGeodesic.geodesic_distance");
173               ErrorLogger.LogException(_className ,"geodesic_coordinate",
174                    new RendererException("Failed inside geodesic_coordinate", exc));
175        }
176        return pt;
177    }
178    /**
179     * Calculates an arc from geodesic point and uses them for the change 1 circular symbols
180     *
181     * @param pPoints array of 3 points, currently the last 2 points are the same. The first point
182     * is the center and the next point defines the radius.
183     *
184     * @return points for the geodesic circle
185     */
186    public static ArrayList<POINT2> GetGeodesicArc(POINT2[] pPoints) {
187        ArrayList<POINT2> pPoints2 = new ArrayList();
188        try {
189            if (pPoints == null) {
190                return null;
191            }
192            if (pPoints.length < 3) {
193                return null;
194            }
195
196            POINT2 ptCenter = new POINT2(pPoints[0]);
197            POINT2 pt1 = new POINT2(pPoints[1]);
198            POINT2 pt2 = new POINT2(pPoints[2]);
199            POINT2 ptTemp = null;
200            ref<double[]> a12b = new ref();
201            double dist2 = 0.0;
202            double dist1 = 0.0;
203            ref<double[]> a12 = new ref();
204            ref<double[]> a21 = new ref();
205            //distance and azimuth from the center to the 1st point
206            dist1 = geodesic_distance(ptCenter, pt1, a12, a21);
207            double saveAzimuth = a21.value[0];
208            //distance and azimuth from the center to the 2nd point
209            dist2 = geodesic_distance(ptCenter, pt2, a12b, a21);
210            //if the points are nearly the same we want 360 degree range fan
211            if (Math.abs(a21.value[0] - saveAzimuth) <= 1) {
212                if (a12.value[0] < 360) {
213                    a12.value[0] += 360;
214                }
215
216                a12b.value[0] = a12.value[0] + 360;
217            }
218
219            ref<double[]> a12c = new ref();
220            int j = 0;
221            if (a12b.value[0] < 0) {
222                a12b.value[0] = 360 + a12b.value[0];
223            }
224            if (a12.value[0] < 0) {
225                a12.value[0] = 360 + a12.value[0];
226            }
227            if (a12b.value[0] < a12.value[0]) {
228                a12b.value[0] = a12b.value[0] + 360;
229            }
230            a12c.value=new double[1];
231            for (j = 0; j <= 100; j++) {
232
233                a12c.value[0] = a12.value[0] + ((double) j / 100.0) * (a12b.value[0] - a12.value[0]);
234                ptTemp = geodesic_coordinate(ptCenter, dist1, a12c.value[0]);
235                pPoints2.add(ptTemp);
236            }
237
238            //if the points are nearly the same we want 360 degree range fan
239            //with no line from the center
240            if (Math.abs(a21.value[0] - saveAzimuth) > 1) {
241                pPoints2.add(ptCenter);
242            }
243
244            if (a12.value[0] < a12b.value[0]) {
245                pPoints2.add(pt1);
246            } else {
247                pPoints2.add(pt2);
248            }
249        } catch (Exception exc) {
250            //clsUtility.WriteFile("Error in mdlGeodesic.GetGeodesicArc");
251               ErrorLogger.LogException(_className ,"GetGeodesicArc",
252                    new RendererException("Failed inside GetGeodesicArc", exc));
253        }
254        return pPoints2;
255    }
256    /**
257     * Calculates the sector points for a sector range fan.
258     *
259     * @param pPoints array of 3 points. The first point
260     * is the center and the next two points define either side of the sector
261     * @param pPoints2 OUT - the calculated geodesic sector points
262     *
263     * @return true if the sector is a circle
264     */
265    public static boolean GetGeodesicArc2(ArrayList<POINT2> pPoints,
266            ArrayList<POINT2> pPoints2) {
267        boolean circle = false;
268        try {
269            POINT2 ptCenter = new POINT2(pPoints.get(0)), pt1 = new POINT2(pPoints.get(1)), pt2 = new POINT2(pPoints.get(2));
270
271            ref<double[]> a12b = new ref();
272            //double dist2 = 0d;
273            double dist1 = 0d;
274            ref<double[]> a12 = new ref();
275            ref<double[]> a21 = new ref();
276            //double lat2c = 0.0;
277            //distance and azimuth from the center to the 1st point
278            //geodesic_distance(lonCenter, latCenter, lon1, lat1, ref dist1, ref a12, ref a21);
279            dist1 = geodesic_distance(ptCenter, pt1, a12, a21);
280            double saveAzimuth = a21.value[0];
281            //distance and azimuth from the center to the 2nd point
282            //geodesic_distance(lonCenter, latCenter, lon2, lat2, ref dist2, ref a12b, ref a21);
283            double dist2 = geodesic_distance(ptCenter, pt2, a12b, a21);
284            //if the points are nearly the same we want 360 degree range fan
285            if (Math.abs(a21.value[0] - saveAzimuth) <= 1) {
286                if (a12.value[0] < 360) {
287                    a12.value[0] += 360;
288                }
289                a12b.value[0] = a12.value[0] + 360;
290                circle = true;
291            }
292
293            //assume caller has set pPoints2 as new Array
294
295            ref<double[]> a12c = new ref();
296            a12c.value = new double[1];
297            int j = 0;
298            POINT2 pPoint = new POINT2();
299            if (a12b.value[0] < 0) {
300                a12b.value[0] = 360 + a12b.value[0];
301            }
302            if (a12.value[0] < 0) {
303                a12.value[0] = 360 + a12.value[0];
304            }
305            if (a12b.value[0] < a12.value[0]) {
306                a12b.value[0] = a12b.value[0] + 360;
307            }
308            for (j = 0; j <= 100; j++) {
309
310                a12c.value[0] = a12.value[0] + ((double) j / 100) * (a12b.value[0] - a12.value[0]);
311                pPoint = geodesic_coordinate(ptCenter, dist1, a12c.value[0]);
312                pPoints2.add(pPoint);
313            }
314        }
315        catch (Exception exc) {
316            //System.out.println(e.getMessage());
317            //clsUtility.WriteFile("Error in mdlGeodesic.GetGeodesicArc2");
318               ErrorLogger.LogException(_className ,"GetGeodesicArc2",
319                    new RendererException("Failed inside GetGeodesicArc2", exc));
320        }
321        return circle;
322    }
323    /**
324     * @deprecated 
325     * returns intersection of two lines, each defined by a point and a bearing
326     * <a rel="license" href="http://creativecommons.org/licenses/by/3.0/"><img alt="Creative Commons License" style="border-width:0" src="http://i.creativecommons.org/l/by/3.0/88x31.png" /></a><br />This work is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by/3.0/">Creative Commons Attribution 3.0 Unported License</a>.
327     * @param p1 1st point
328     * @param brng1 first line bearing in degrees from true north
329     * @param p2 2nd point
330     * @param brng2 2nd point bearing in degrees from true north
331     * @return
332     */
333    public static POINT2 IntersectLines(POINT2 p1,
334            double brng1, 
335            POINT2 p2,
336            double brng2) {
337        POINT2 ptResult = null;
338        try {
339            double lat1 = DegToRad(p1.y);//p1._lat.toRad();
340            double lon1 = DegToRad(p1.x);//p1._lon.toRad();
341            double lat2 = DegToRad(p2.y);//p2._lat.toRad();
342            double lon2 = DegToRad(p2.x);//p2._lon.toRad();
343            double brng13 = DegToRad(brng1);//brng1.toRad();
344            double brng23 = DegToRad(brng2);//brng2.toRad();
345            double dLat = lat2 - lat1;
346            double dLon = lon2 - lon1;
347
348
349            double dist12 = 2 * Math.asin(Math.sqrt(Math.sin(dLat / 2) * Math.sin(dLat / 2) +
350                    Math.cos(lat1) * Math.cos(lat2) * Math.sin(dLon / 2) * Math.sin(dLon / 2)));
351
352            if (dist12 == 0) {
353                return null;
354            }
355
356            double brngA = Math.acos((Math.sin(lat2) - Math.sin(lat1) * Math.cos(dist12)) /
357                    (Math.sin(dist12) * Math.cos(lat1)));
358
359            if (Double.isNaN(brngA)) {
360                brngA = 0;  // protect against rounding
361            }
362            double brngB = Math.acos((Math.sin(lat1) - Math.sin(lat2) * Math.cos(dist12)) /
363                    (Math.sin(dist12) * Math.cos(lat2)));
364
365            double brng12 = 0, brng21 = 0;
366            if (Math.sin(lon2 - lon1) > 0) {
367                brng12 = brngA;
368                brng21 = 2 * Math.PI - brngB;
369            } else {
370                brng12 = 2 * Math.PI - brngA;
371                brng21 = brngB;
372            }
373
374            double alpha1 = (brng13 - brng12 + Math.PI) % (2 * Math.PI) - Math.PI;  // angle 2-1-3
375            double alpha2 = (brng21 - brng23 + Math.PI) % (2 * Math.PI) - Math.PI;  // angle 1-2-3
376
377            if (Math.sin(alpha1) == 0 && Math.sin(alpha2) == 0) {
378                return null;  // infinite intersections
379            }
380            if (Math.sin(alpha1) * Math.sin(alpha2) < 0) {
381                return null;       // ambiguous intersection
382            }
383            //alpha1 = Math.abs(alpha1);
384            //alpha2 = Math.abs(alpha2);  // ... Ed Williams takes abs of alpha1/alpha2, but seems to break calculation?
385            double alpha3 = Math.acos(-Math.cos(alpha1) * Math.cos(alpha2) +
386                    Math.sin(alpha1) * Math.sin(alpha2) * Math.cos(dist12));
387
388            double dist13 = Math.atan2(Math.sin(dist12) * Math.sin(alpha1) * Math.sin(alpha2),
389                    Math.cos(alpha2) + Math.cos(alpha1) * Math.cos(alpha3));
390
391            double lat3 = Math.asin(Math.sin(lat1) * Math.cos(dist13) +
392                    Math.cos(lat1) * Math.sin(dist13) * Math.cos(brng13));
393            double dLon13 = Math.atan2(Math.sin(brng13) * Math.sin(dist13) * Math.cos(lat1),
394                    Math.cos(dist13) - Math.sin(lat1) * Math.sin(lat3));
395            double lon3 = lon1 + dLon13;
396            lon3 = (lon3 + Math.PI) % (2 * Math.PI) - Math.PI;  // normalise to -180..180º
397
398            //return new POINT2(lat3.toDeg(), lon3.toDeg());
399            ptResult = new POINT2(RadToDeg(lon3), RadToDeg(lat3));
400
401        } catch (Exception exc) {
402            ErrorLogger.LogException(_className, "IntersectLines",
403                    new RendererException("Failed inside IntersectLines", exc));
404        }
405        return ptResult;
406    }
407    /**
408     * Normalizes geo points for arrays which span the IDL
409     *
410     * @param geoPoints
411     * @return
412     */
413    public static ArrayList<POINT2> normalize_points(ArrayList<POINT2> geoPoints) {
414        ArrayList<POINT2> normalizedPts = null;
415        try {
416            if (geoPoints == null || geoPoints.isEmpty()) {
417                return normalizedPts;
418            }
419
420            int j = 0;
421            double minx = geoPoints.get(0).x;
422            double maxx = minx;
423            boolean spansIDL = false;
424            POINT2 pt = null;
425            int n=geoPoints.size();
426            //for (j = 1; j < geoPoints.size(); j++) 
427            for (j = 1; j < n; j++) 
428            {
429                pt = geoPoints.get(j);
430                if (pt.x < minx) {
431                    minx = pt.x;
432                }
433                if (pt.x > maxx) {
434                    maxx = pt.x;
435                }
436            }
437            if (maxx - minx > 180) {
438                spansIDL = true;
439            }
440
441            if (!spansIDL) {
442                return geoPoints;
443            }
444
445            normalizedPts = new ArrayList();
446            n=geoPoints.size();
447            //for (j = 0; j < geoPoints.size(); j++) 
448            for (j = 0; j < n; j++) 
449            {
450                pt = geoPoints.get(j);
451                if (pt.x < 0) {
452                    pt.x += 360;
453                }
454                normalizedPts.add(pt);
455            }
456        } catch (Exception exc) {
457            ErrorLogger.LogException(_className, "normalize_pts",
458                    new RendererException("Failed inside normalize_pts", exc));
459        }
460        return normalizedPts;
461    }
462
463    /**
464     * calculates the geodesic MBR, intended for regular shaped areas
465     *
466     * @param geoPoints
467     * @return
468     */
469    public static Rectangle2D.Double geodesic_mbr(ArrayList<POINT2> geoPoints) {
470        Rectangle2D.Double rect2d = null;
471        try {
472            if (geoPoints == null || geoPoints.isEmpty()) {
473                return rect2d;
474            }
475            
476            ArrayList<POINT2>normalizedPts=normalize_points(geoPoints);
477            double ulx=normalizedPts.get(0).x;
478            double lrx=ulx;
479            double uly=normalizedPts.get(0).y;
480            double lry=uly;
481            int j=0;
482            POINT2 pt=null;
483            int n=normalizedPts.size();
484            //for(j=1;j<normalizedPts.size();j++)
485            for(j=1;j<n;j++)
486            {
487                pt=normalizedPts.get(j);
488                if(pt.x<ulx)
489                    ulx=pt.x;
490                if(pt.x>lrx)
491                    lrx=pt.x;
492            
493                if(pt.y>uly)
494                    uly=pt.y;
495                if(pt.y<lry)
496                    lry=pt.y;
497            }
498            POINT2 ul=new POINT2(ulx,uly);
499            POINT2 ur=new POINT2(lrx,uly);
500            POINT2 lr=new POINT2(lrx,lry);
501            double width=geodesic_distance(ul,ur,null,null);
502            double height=geodesic_distance(ur,lr,null,null);
503            rect2d=new Rectangle2D.Double(ulx,uly,width,height);
504        } catch (Exception exc) {
505            ErrorLogger.LogException(_className, "geodesic_mbr",
506                    new RendererException("Failed inside geodesic_mbr", exc));
507        }
508        return rect2d;
509    }
510
511    /**
512     * Currently used by AddModifiers for greater accuracy on center labels
513     *
514     * @param geoPoints
515     * @return
516     */
517    public static POINT2 geodesic_center(ArrayList<POINT2> geoPoints) {
518        POINT2 pt = null;
519        try {
520            if(geoPoints==null || geoPoints.isEmpty())
521                return pt;
522            
523            Rectangle2D.Double rect2d=geodesic_mbr(geoPoints);
524            double deltax=rect2d.getWidth()/2;
525            double deltay=rect2d.getHeight()/2;
526            POINT2 ul=new POINT2(rect2d.x,rect2d.y);
527            //first walk east by deltax
528            POINT2 ptEast=geodesic_coordinate(ul,deltax,90);
529            //next walk south by deltay;
530            pt=geodesic_coordinate(ptEast,deltay,180);
531            
532        } catch (Exception exc) {
533            ErrorLogger.LogException(_className, "geodesic_center",
534                    new RendererException("Failed inside geodesic_center", exc));
535        }
536        return pt;
537    }
538    /**
539     * rotates a point from a center point in degrees
540     * @param ptCenter center point to rotate about
541     * @param ptRotate point to rotate
542     * @param rotation rotation angle in degrees
543     * @return 
544     */
545    private static POINT2 geoRotatePoint(POINT2 ptCenter, POINT2 ptRotate, double rotation)
546    {
547        try
548        {
549            double bearing=GetAzimuth(ptCenter, ptRotate);
550            double dist=geodesic_distance(ptCenter,ptRotate,null,null);
551            return geodesic_coordinate(ptCenter,dist,bearing+rotation);
552        }
553        catch (Exception exc) {
554            ErrorLogger.LogException(_className, "geoRotatePoint",
555                    new RendererException("Failed inside geoRotatePoint", exc));
556        }
557        return null;
558    }
559    /**
560     * Calculates points for a geodesic ellipse and rotates the points by rotation
561     * @param ptCenter
562     * @param majorRadius
563     * @param minorRadius
564     * @param rotation  rotation angle in degrees
565     * @return 
566     */
567    public static POINT2[] getGeoEllipse(POINT2 ptCenter, double majorRadius, double minorRadius, double rotation)
568    {        
569        POINT2[]pEllipsePoints=null;
570        try
571        {
572            pEllipsePoints=new POINT2[37];
573            //int l=0;
574            POINT2 pt=null;            
575            double dFactor, azimuth=0,a=0,b=0,dist=0,bearing=0;
576            POINT2 ptLongitude=null,ptLatitude=null;
577            for (int l = 1; l < 37; l++)
578            {
579                dFactor = (10.0 * l) * Math.PI / 180.0;                
580                a=majorRadius * Math.cos(dFactor);
581                b=minorRadius * Math.sin(dFactor);
582                //dist=Math.sqrt(a*a+b*b);
583                //azimuth = (10.0 * l);// * Math.PI / 180.0;  
584                //azimuth=90-azimuth;
585                //pt = geodesic_coordinate(ptCenter,dist,azimuth);                
586                //pt = geodesic_coordinate(ptCenter,dist,azimuth);                
587                ptLongitude=geodesic_coordinate(ptCenter,a,90);
588                ptLatitude=geodesic_coordinate(ptCenter,b,0);
589                //pt=new POINT2(ptLatitude.x,ptLongitude.y);
590                pt=new POINT2(ptLongitude.x,ptLatitude.y);
591                //pEllipsePoints[l-1]=pt;
592                pEllipsePoints[l-1]=geoRotatePoint(ptCenter,pt,-rotation);
593            }            
594            pEllipsePoints[36]=new POINT2(pEllipsePoints[0]);
595        }
596        catch(Exception exc)
597        {
598            ErrorLogger.LogException(_className, "GetGeoEllipse",
599                    new RendererException("GetGeoEllipse", exc));
600        }
601        return pEllipsePoints;
602    }
603}