transforms.cc

   1 /*
   2  * Transformations
   3  */
   4
   5 #include <math.h>
   6 #include <stdio.h>
   7
   8 #include "transforms.hh"
   9
  10 RotationTransform::RotationTransform(double theta) {
  11   s= sin(theta); c= cos(theta);
  12 }
  13
  14 Point XZRotationTransform::operator()(Point in) {
  15   return Point(in[0]*c + in[2]*s,
  16                in[1],
  17                in[2]*c - in[0]*s);
  18 }
  19
  20 Point YZRotationTransform::operator()(Point in) {
  21   return Point(in[0],
  22                in[1]*c + in[2]*s,
  23                in[2]*c - in[1]*s);
  24 }
  25
  26 TearOpenTransform::TearOpenTransform(double amt, double apy) {
  27   factor= pow(10.0,amt);
  28   apexy= apy;
  29 }
  30
  31 Point TearOpenTransform::operator()(Point in) {
  32   double w= apexy - in[1];
  33   double r= sqrt(in[2]*in[2] + w*w);
  34   double theta= r>0.0 ? factor*atan2(in[2],w) : 0.0;
  35   double ny= apexy - r*cos(theta);
  36   double nz= r*sin(theta);
  37   return Point(in[0],ny,nz);
  38 }
  39
  40 BulgeTransform::BulgeTransform(double amt, double stddev, double apy) {
  41   amount=amt; apexy=apy;
  42   variance= stddev*stddev;
  43 }
  44
  45 Point BulgeTransform::operator()(Point in) {
  46   double w= in[1] - apexy;
  47   double r2= w*w + in[2]*in[2];
  48   double nror= 1.0 + amount*exp(-r2/variance);
  49   double ny= w*nror + apexy;
  50   return Point(in[0],ny,in[2]);
  51 }
  52
  53 Point ShearTransform::operator()(Point in) {
  54   double dx= (in[0] - epicentre[0]) / deviation[0];
  55   double dy= (in[1] - epicentre[1]) / deviation[1];
  56   double dz= (in[2] - epicentre[2]) / deviation[2];
  57   double r2= dx*dx + dy*dy + dz*dz;
  58   double amount= exp(-r2);
  59   return in + magnitude*amount;
  60 }
  61
  62 Point TwistTransform::operator()(Point in) {
  63   double w= in[1] - apexy;
  64   double theta= factor*in[2];
  65   double c= cos(theta), s= sin(theta);
  66   double nx= in[0]*c - w*s;
  67   double ny= w*c + in[0]*s + apexy;
  68   return Point(nx,ny,in[2]);
  69 }
  70
  71 Point ScaleTransform::operator()(Point in) {
  72   return in*scale;
  73 }
  74
  75 Point MagicTransform::operator()(Point in) {
  76   double nz= in[2]*zfactor;
  77   double r= sqrt(in[0]*in[0] + in[1]*in[1]);
  78   double nx= xfactor * (r - xkonstant);
  79   double nr= yoffset + yfactor*in[1];
  80   double ny2=  nr*nr- nz*nz;
  81 //fprintf(stderr,"M %7.4f,%7.4f,%7.4f r=%7.4f nx=%7.4f nz=%7.4f nr=%7.4f ny2=%7.4f\n",
  82 //        in[0],in[1],in[2], r,nx,nz,nr,ny2);
  83   double ny= sqrt(ny2);
  84   return Point(nx,ny,nz);
  85 }