trying great circles: seems to be an error

1 files changed, 110 insertions, 55 deletions

diff --git a/vodaviz.pde b/vodaviz.pde
index 9a4272f..f29001c 100644
--- a/vodaviz.pde
+++ b/vodaviz.pde
@@ -42,27 +42,66 @@
  O- The orientation of the camera (represented, for instance, by Tait–Bryan angles).
  e- the viewer's position relative to the display surface.[1]
 
+
+http://williams.best.vwh.net/avform.htm#Crs
+==good stuff
+
+Intermediate points on a great circle
+Here we find points (lat,lon) a given fraction of the distance (d) between them.
+Suppose the starting point is (lat1,lon1) and the final point (lat2,lon2) and we want the point
+a fraction f along the great circle route. f=0 is point 1. f=1 is point 2. The two points cannot
+be antipodal ( i.e. lat1+lat2=0 and abs(lon1-lon2)=pi) because then the route is undefined.
+The intermediate latitude and longitude is then given by:
+
+        A=sin((1-f)*d)/sin(d)
+        B=sin(f*d)/sin(d)
+        x = A*cos(lat1)*cos(lon1) +  B*cos(lat2)*cos(lon2)
+        y = A*cos(lat1)*sin(lon1) +  B*cos(lat2)*sin(lon2)
+        z = A*sin(lat1)           +  B*sin(lat2)
+        lat=atan2(z,sqrt(x^2+y^2))
+        lon=atan2(y,x)
+        
+        //find 
+        
+        RPoint sp,ep,mp;
+        float f=0.8; //fraction along path;
+        float d=acos(sin(sp.y)*sin(ep.y)+cos(sp.y)*cos(ep.y)*cos(sp.x-ep.x));
+        float A=sin((1-f)*d)/sin(d);
+        float B=sin(f*d)/sin(d);
+        float x = A*cos(lat1)*cos(lon1) +  B*cos(lat2)*cos(lon2);
+        float y = A*cos(lat1)*sin(lon1) +  B*cos(lat2)*sin(lon2);
+        float z = A*sin(lat1)           +  B*sin(lat2);
+        mp.y=atan2(z,sqrt(x^2+y^2));
+        mp.x=atan2(y,x);
+
  */
 
 import processing.pdf.*;
 import geomerative.*;
 
-class pointMapper {
+class pointNormalise {
+  //take pixel coords and turn into lat/lng radians
   float xo,xs,yo,ys;
-  pointMapper(float _xo,float _xs,float _yo,float _ys) {
+  pointNormalise(float _xo,float _xs,float _yo,float _ys) {
     xo=_xo;
     xs=_xs;
     yo=_yo;
     ys=_ys;
   }
-  RPoint map(float _x,float _y) {
+  RPoint alise(float _x,float _y) {
         //formalise
-        float px=((_x-xo)/xs); //26)/736);
-        float py=((_y-yo)/ys); //90)/390);
+        float px=((_x-xo)/xs)*PI*2; //26)/736);
+        float py=((_y-yo)/ys)*PI; //90)/390);
+        
+        return new RPoint(px,py);
+  }
+}
+
+class pointMapper {
         
         //transform
-        px=(((px-0.5)*0.25)+.25); //front half of sphere
-        py=((py-0.5)*0.3)+0.5;
+        float px=(((_x-0.5)*0.25)+.25); //front half of sphere
+        float py=((_y-0.5)*0.3)+0.5;
         
         //map to 3D curve
         float r=getHeight()*4;
@@ -82,15 +121,16 @@ class pointMapper {
 }
 
 RShape shp;
+pointNormalise pnorm;
 pointMapper shpmap;
- //RRectangle cb;
-int whichChild,numChildren;
 
 String mode;
 
 csvloader data;
 calldata calls;
 
+country Ireland;
+
 void setup(){
   println("vodaviz v0.11");
   RG.init(this);
@@ -103,69 +143,84 @@ void setup(){
   
   shp = RG.loadShape("countries_named_mercator.svg"); //test_drawing.svg"); //world_countries_outlines_split.svg");
   
-  shpmap = new pointMapper(26,736,90,390); 
+  pnorm = new pointNormalise(26,736,90,390); 
+  shpmap = new pointMapper(); 
+  
   RG.ignoreStyles();
   println("loaded svg in "+((millis()-m)*.001)+" seconds");
-  numChildren=shp.children.length;
-  whichChild=0;
-
+  Ireland=new country("Ireland",0,0,shp.children[0]);
   data=new csvloader("calls.csv");
-  calls=new calldata(data.data);
+  calls=new calldata(data.data,shp);
  
-}
-
-void draw() {
   background(0,0,0);
   noFill();
   stroke(255);
   strokeWeight(.02);
-  //RG.shape(shp);
-  whichChild=(whichChild+1)%numChildren;
-  //println("child "+whichChild+" of "+numChildren);
-  
-  //cb=shp.children[whichChild].getBounds();
   
-  float xWrap=1.0;
-  
-  //plot shapes at points
-  for (int i=0;i<numChildren;i++) { 
-    String cmd="";
-    for (int j=0;j<shp.children[i].paths.length;j++) { //
-      //shp.children[i].paths[j].draw();
-      
-      //public boolean contains(RPoint p)
-      
-      
-      /*
-      //find a random point inside shape
-      RPoint tl=shp.children[i].getTopLeft();
-      RPoint br=shp.children[i].getBottomRight();
-      RPoint np;
-      boolean found =false;
-      while (!found) {
-        np=new RPoint(random(tl.x,br.x),random(tl.y,br.y));
-        found = shp.children[i].contains(np); //paths[j].
-        println(shp.children[i].name+"; trying "+np.x+","+np.y+" between "+tl.x+","+tl.y+" and "+br.x+","+br.y+" :"+found);
-      }
-      */
-      
-      
-      
+  //pick a random colour
+  stroke(random(200)+50, random(200)+50,random(200)+50);
+  println("plotting "+calls.countries.get(0).name+": "+calls.countries.get(0).calls+" calls");
+}
+
+int i=0;
+int j=0;
+
+void draw() {
+ 
+  //for (int i=0;i<calls.countries.size();i++) { 
+
+       /*
       beginShape();
-      for (int k=0;k<shp.children[i].paths[j].commands.length;k++) {
-        RPoint sp=shp.children[i].paths[j].commands[k].startPoint;
+      for (int k=0;k<calls.getcountry(i).outline.paths[0].commands.length;k++) {
+        RPoint sp=calls.getcountry(i).outline.paths[0].commands[k].startPoint;
         RPoint dp=shpmap.map(sp.x,sp.y);
         
         vertex(dp.x+(getWidth()/2),dp.y+(getHeight()/2)); //,z); //z);
-        //if (i==0&&j==0) cmd+=(x+(getWidth()/2))+","+(y+(getHeight()/2))+","+z+" ";
       }
       endShape();
+      */
       
-    }
-    if (cmd!="") println(cmd);
+      //pick a random colour
+      //stroke(random(100)+100, random(100)+100,random(100)+100);
+      //for (int j=0;j<calls.countries.get(i).calls;j++) {
+        RPoint sp=pnorm.alise(calls.countries.get(i).getpoint(100));
+        RPoint ep=pnorm.alise(Ireland.getpoint(100));
+        if (sp.x>0&&ep.x>0){
+        RPoint Sp=shpmap.map(sp.x,sp.y);
+        RPoint Ep=shpmap.map(ep.x,ep.y);
+          
+        //find point along path
+        
+        float f=0.8; //fraction along path;
+        float d=acos(sin(sp.y)*sin(ep.y)+cos(sp.y)*cos(ep.y)*cos(sp.x-ep.x));
+        float A=sin((1-f)*d)/sin(d);
+        float B=sin(f*d)/sin(d);
+        float x = A*cos(sp.y)*cos(sp.x) +  B*cos(ep.y)*cos(ep.x);
+        float y = A*cos(sp.y)*sin(sp.x) +  B*cos(ep.y)*sin(ep.x);
+        float z = A*sin(sp.y)           +  B*sin(ep.y);
+        RPoint mp=new RPoint(atan2(y,x),atan2(z,sqrt(pow(x,2)+pow(y,2))));
+        
+        RPoint Mp=shpmap.map(mp.x,mp.y);
+          beginShape();
+            line(Sp.x+(getWidth()/2),Sp.y+(getHeight()/2),Mp.x+(getWidth()/2),Mp.y+(getHeight()/2));
+            line(Mp.x+(getWidth()/2),Mp.y+(getHeight()/2),Ep.x+(getWidth()/2),Ep.y+(getHeight()/2));
+          endShape();
+        }
+      //}
+      j=j+1;
+      if (j>(calls.countries.get(i).calls*.01)){
+        i++;
+        stroke(random(200)+50, random(200)+50,random(200)+50);
+        println("plotting "+calls.countries.get(i).name+": "+calls.countries.get(i).calls+" calls");
+        j=0;
+      }
+
+  //}
+  if (i==calls.countries.size()-1) {
+    println("finished");
+   noLoop();
   }
- noLoop();
- if (mode=="PDF") exit();
+ //if (mode=="PDF") exit();
 }
 
 void mousePressed() {