Repatch (from linux) of OSX IRAF

1 files changed, 214 insertions, 0 deletions

diff --git a/noao/rv/complex.x b/noao/rv/complex.x
new file mode 100644
index 00000000..68fb06e2
--- /dev/null
+++ b/noao/rv/complex.x
@@ -0,0 +1,214 @@
+include <math.h>
+include <mach.h>
+
+# COMPLEX.X  - File containing utility routines for complex arithmetic.
+
+# CX_ADD - Addition of complex numbers.
+
+procedure cx_add (ar, ai, br, bi, cr, ci)
+
+real	ar, ai					#I First number
+real	br, bi					#I Second number
+real	cr, ci					#O Computed value
+
+begin
+	cr = ar + br
+	ci = ai + bi
+end
+
+
+# CX_SUB - Subtraction of complex numbers.
+
+procedure cx_sub (ar, ai, br, bi, cr, ci)
+
+real	ar,ai					#I First number
+real	br,bi					#I Second number
+real	cr,ci					#O Computed value
+
+begin
+	cr = ar - br
+	ci = ai - bi
+end
+
+
+# CX_MUL - Multiplication of complex numbers.
+
+procedure cx_mul (ar, ai, br, bi, cr, ci)
+
+real	ar,ai					#I First number
+real	br,bi					#I Second number
+real	cr,ci					#O Computed value
+
+begin
+	cr = ar*br - ai*bi
+	ci = ai*br + ar*bi
+end
+
+
+# CX_DIV - Division of complex numbers.
+
+procedure cx_div  (ar, ai, br, bi, cr, ci)
+
+real	ar,ai					#I First number
+real	br,bi					#I Second number
+real	cr,ci					#O Computed value
+
+real	r, den
+
+begin
+	if (br == 0.0 && bi == 0.0) {		# Trap divide by zero
+	    cr = 0.0
+	    ci = 0.0
+	    return
+	}
+
+	if (abs(br) >= abs(bi)) {
+	    r = bi / br
+	    den = br + r*bi
+	    cr = (ar + r*ai) / den
+	    ci = (ai - r*ar) / den
+	} else {
+	    r = br / bi
+	    den = bi + r*br
+	    cr = (ar*r + ai) / den
+	    ci = (ai*r - ar) / den
+	}
+end
+
+
+# CX_ABS - Absolute value of complex numbers.
+
+real procedure cx_abs (ar, ai)
+
+real	ar, ai					#I First number
+
+real	x, y, ans, temp
+
+begin
+	x = abs (ar)
+	y = abs (ai)
+	if (x == 0.0)
+	    ans = y
+	else if (y == 0.0)
+	    ans = x
+	else if (x > y) {
+	    temp = y / x
+	    ans = x * sqrt (1.0 + temp*temp)
+	} else {
+	    temp = x / y
+	    ans = y * sqrt (1.0 + temp*temp)
+	}
+
+	return (ans)
+end
+
+
+# CX_CONJG - Complex conjugate.
+
+procedure cx_conjg (ar, ai, br, bi)
+
+real	ar,ai					#I First number
+real	br,bi					#I Conjugate
+
+begin
+	br =   ar
+	bi = - ai
+end
+
+
+# CX_SQRT - Square root of complex numbers.
+
+procedure cx_sqrt (ar, ai, br, bi)
+
+real	ar, ai					#I First number
+real	br, bi					#I Square root
+
+real	x, y, w, r
+
+begin
+	if (ar == 0.0 && ai == 0.0) {
+	    br = 0.0
+	    bi = 0.0
+	} else {
+	    x = abs (ar)
+	    y = abs (ai)
+	    if (x >= y) {
+		r = y / x
+		w = sqrt (x) * sqrt (0.5*(1.0+sqrt(1.0+r*r)))
+	    } else {
+		r = x / y
+		w = sqrt (y) * sqrt (0.5*(1.0+sqrt(1.0+r*r)))
+	    }
+	    if (ar >= 0.0) {
+		br = w
+		bi = ai / (2.0*w)
+	    } else {
+		if (ai >= 0)
+		    bi =  w
+		else
+		    bi = -w
+		br = ai / (2.0*bi) 
+	    }
+	}
+end
+
+
+# CEXP1 - Complex exponentiation routine.
+
+procedure cexp1 (a, b, dr, di) 
+
+real	a				#I Real part of argument
+real	b				#I Complex part of argument
+real	dr, di				#O Resultant real/imaginary components
+
+begin
+	if (a > log(MAX_REAL)) {
+	    dr = 0.0
+	    di = 0.0
+	} else
+	    call cx_div (cos(b), sin(b), exp(a), 0.0, dr, di)
+end 
+
+
+# CX_PAK - Pack two real arrays of an FFT into one real FFT array.
+# The array `fft' must be dimensioned to at least 2*fnpts elements.
+
+procedure cx_pak (creal, cimg, fft, fnpts) 
+
+real	creal[fnpts], cimg[fnpts]	#I Real/Img complex components
+real	fft[ARB]			#O Output 'real' array
+int	fnpts				#I Npts in array
+
+int	i,j
+
+begin
+	j = 1
+	do i = 1, fnpts {
+	    fft[j] = creal[i]
+	    j = j + 1
+	    fft[j] = cimg[i]
+	    j = j + 1
+	}
+end
+
+
+# CX_UNPAK - Unpack one real FFT array into two component real arrays.
+# The array `fft' must be dimensioned to at least 2*fnpts elements.
+
+procedure cx_unpak (fft, creal, cimg, fnpts) 
+
+real	fft[ARB]			#O Output 'real' array
+real	creal[fnpts], cimg[fnpts]	#I Real/Img complex components
+int	fnpts				#I Npts in array
+
+int	i,j
+
+begin
+	j = 1
+	do i = 1, fnpts {
+	    creal[i] = fft[j]
+	    j = j + 1
+	    cimg[i] = fft[j]
+	    j = j + 1
+	}
+end