gonum.org/v1/gonum@v0.14.0/num/dualcmplx/dual_example_test.go (about)

     1  // Copyright ©2018 The Gonum Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package dualcmplx_test
     6  
     7  import (
     8  	"fmt"
     9  	"math"
    10  
    11  	"gonum.org/v1/gonum/floats/scalar"
    12  	"gonum.org/v1/gonum/num/dualcmplx"
    13  )
    14  
    15  // point is a 2-dimensional point/vector.
    16  type point struct {
    17  	x, y float64
    18  }
    19  
    20  // raise raises the dimensionality of a point to a complex.
    21  func raise(p point) complex128 {
    22  	return complex(p.x, p.y)
    23  }
    24  
    25  // raiseDual raises the dimensionality of a point to a dual complex number.
    26  func raiseDual(p point) dualcmplx.Number {
    27  	return dualcmplx.Number{
    28  		Real: 1,
    29  		Dual: complex(p.x, p.y),
    30  	}
    31  }
    32  
    33  // transform performs the transformation of p by the given dual complex numbers.
    34  // The transformations are normalized to unit vectors.
    35  func transform(p point, by ...dualcmplx.Number) point {
    36  	if len(by) == 0 {
    37  		return p
    38  	}
    39  
    40  	// Ensure the modulus of by is correctly scaled.
    41  	for i := range by {
    42  		if len := dualcmplx.Abs(by[i]); len != 1 {
    43  			by[i].Real *= complex(1/len, 0)
    44  		}
    45  	}
    46  
    47  	// Perform the transformations.
    48  	z := by[0]
    49  	for _, o := range by[1:] {
    50  		z = dualcmplx.Mul(o, z)
    51  	}
    52  	pp := dualcmplx.Mul(dualcmplx.Mul(z, raiseDual(p)), dualcmplx.Conj(z))
    53  
    54  	// Extract the point.
    55  	return point{x: real(pp.Dual), y: imag(pp.Dual)}
    56  }
    57  
    58  func Example() {
    59  	// Translate a 1×1 square by [3, 4] and rotate it 90° around the
    60  	// origin.
    61  	fmt.Println("square:")
    62  
    63  	// Construct a displacement.
    64  	displace := dualcmplx.Number{
    65  		Real: 1,
    66  		Dual: 0.5 * raise(point{3, 4}),
    67  	}
    68  
    69  	// Construct a rotation.
    70  	alpha := math.Pi / 2
    71  	rotate := dualcmplx.Number{Real: complex(math.Cos(alpha/2), math.Sin(alpha/2))}
    72  
    73  	for i, p := range []point{
    74  		{x: 0, y: 0},
    75  		{x: 0, y: 1},
    76  		{x: 1, y: 0},
    77  		{x: 1, y: 1},
    78  	} {
    79  		pp := transform(p,
    80  			displace, rotate,
    81  		)
    82  
    83  		// Clean up floating point error for clarity.
    84  		pp.x = scalar.Round(pp.x, 2)
    85  		pp.y = scalar.Round(pp.y, 2)
    86  
    87  		fmt.Printf(" %d %+v -> %+v\n", i, p, pp)
    88  	}
    89  
    90  	// Rotate a line segment 90° around its lower end [2, 2].
    91  	fmt.Println("\nline segment:")
    92  
    93  	// Construct a displacement to the origin from the lower end...
    94  	origin := dualcmplx.Number{
    95  		Real: 1,
    96  		Dual: 0.5 * raise(point{-2, -2}),
    97  	}
    98  	// ... and back from the origin to the lower end.
    99  	replace := dualcmplx.Number{
   100  		Real: 1,
   101  		Dual: -origin.Dual,
   102  	}
   103  
   104  	for i, p := range []point{
   105  		{x: 2, y: 2},
   106  		{x: 2, y: 3},
   107  	} {
   108  		pp := transform(p,
   109  			origin,  // Displace to origin.
   110  			rotate,  // Rotate around axis.
   111  			replace, // Displace back to original location.
   112  		)
   113  
   114  		// Clean up floating point error for clarity.
   115  		pp.x = scalar.Round(pp.x, 2)
   116  		pp.y = scalar.Round(pp.y, 2)
   117  
   118  		fmt.Printf(" %d %+v -> %+v\n", i, p, pp)
   119  	}
   120  
   121  	// Output:
   122  	//
   123  	// square:
   124  	//  0 {x:0 y:0} -> {x:-4 y:3}
   125  	//  1 {x:0 y:1} -> {x:-5 y:3}
   126  	//  2 {x:1 y:0} -> {x:-4 y:4}
   127  	//  3 {x:1 y:1} -> {x:-5 y:4}
   128  	//
   129  	// line segment:
   130  	//  0 {x:2 y:2} -> {x:2 y:2}
   131  	//  1 {x:2 y:3} -> {x:1 y:2}
   132  }