chiark / gitweb /
found in ijackson@chiark:things/moebius.old-before-cvs
[moebius.git] / 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 MagicTransform::operator()(Point in) {
72   double nz= in[2]*zfactor;
73   double r= sqrt(in[0]*in[0] + in[1]*in[1]);
74   double nx= xfactor * (r - xkonstant);
75   double nr= yoffset + yfactor*in[1];
76   double ny2=  nr*nr- nz*nz;
77 //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",
78 //        in[0],in[1],in[2], r,nx,nz,nr,ny2);
79   double ny= sqrt(ny2);
80   return Point(nx,ny,nz);
81 }