//vodafone call visualisation
//
//Tim Redfern Jan 2012
//
//scale of wall 10.57x2.8m - some parts are obscured
//
//416 x 110 in. ~ 62400x16500 ~ 1029 mpx ~ 2945 MB ram to open
//
//CousinMarriageWorld.svg, nominally 940 × 470 pixels, file size: 1.99 MB
//World_map_(Miller_cylindrical_projection,_blank).svg‎, nominally 634 × 477 pixels, file size: 1.84 MB
//wikimedia
//
//http://www.vectortemplates.com/vector-world-map.php

//requirements - 
//must be able to transform points to a new projection
// maybe not through normal shape library
//must be able to find a random point within the shape
// 

//250112 so far so good, can iterate children.
//safest thing is to establish that I can place a point as being inside or outside of a shape
//(are they closed properly!)

//If we can find points in them then we can start the main task:
//connecting parts of the SVG with rows in the spreadsheet

//maybe next first I should look at projections (to keep them happy)
//I have a map with rectangular coords,
//is easy to enough to play with the projection

// public boolean contains(RGeomElem shp)
// public geomerative.RRectangle getBounds()
// --> this is kind of hidden. why? anyway, it should be possible to use it
//another option is to use the OUTLINE of the country (as a purely geometric thing).

/*
 http://en.wikipedia.org/wiki/3D_projection
 
 a- the 3D position of a point A that is to be projected.
 c- the 3D position of a point C representing the camera.
 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 pointNormalise {
  //take pixel coords and turn into lat/lng radians
  float xo,xs,yo,ys;
  pointNormalise(float _xo,float _xs,float _yo,float _ys) {
    xo=_xo;
    xs=_xs;
    yo=_yo;
    ys=_ys;
  }
  RPoint alise(float _x,float _y) {
        //formalise
        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
        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;
        float x=r*cos(px*PI*2)*(sin(py*PI)) *2;
        float z=r*sin(px*PI*2)*(sin(py*PI))-985;
        float y=r*cos(py*PI)+(getHeight()*0.5)-210;  
        
        //camera at 0,0,0 
        //screen plane at 0,0,100
        //(dx-ex)(ez/dz)
        //(dy-ey)(ez/dz)
        
        //rotate camera

        return new RPoint(x*(100/z),y*(100/z));
    }
}

RShape shp;
pointNormalise pnorm;
pointMapper shpmap;

String mode;

csvloader data;
calldata calls;

country Ireland;

void setup(){
  println("vodaviz v0.11");
  RG.init(this);

  //mode="PDF";
  size(832,220); //,PDF, "testoutput.pdf"); //P3D); //832,220); //nb pdf is 800x600
  
  smooth();
  float m = millis();
  
  shp = RG.loadShape("countries_named_mercator.svg"); //test_drawing.svg"); //world_countries_outlines_split.svg");
  
  pnorm = new pointNormalise(26,736,90,390); 
  shpmap = new pointMapper(); 
  
  RG.ignoreStyles();
  println("loaded svg in "+((millis()-m)*.001)+" seconds");
  Ireland=new country("Ireland",0,0,shp.children[0]);
  data=new csvloader("calls.csv");
  calls=new calldata(data.data,shp);
 
  background(0,0,0);
  noFill();
  stroke(255);
  strokeWeight(.02);
  
  //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<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);
      }
      endShape();
      */
      
      //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();
  }
 //if (mode=="PDF") exit();
}

void mousePressed() {
  
}