8 #include "transforms.hh"
10 RotationTransform::RotationTransform(double theta) {
11 s= sin(theta); c= cos(theta);
14 Point XZRotationTransform::operator()(Point in) {
15 return Point(in[0]*c + in[2]*s,
20 Point YZRotationTransform::operator()(Point in) {
26 TearOpenTransform::TearOpenTransform(double amt, double apy) {
27 factor= pow(10.0,amt);
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);
40 BulgeTransform::BulgeTransform(double amt, double stddev, double apy) {
41 amount=amt; apexy=apy;
42 variance= stddev*stddev;
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]);
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;
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]);
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);
80 return Point(nx,ny,nz);