forked from foxsi/calsoft
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcmreplicate.pro
executable file
·146 lines (132 loc) · 4.88 KB
/
cmreplicate.pro
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
;+
; NAME:
; CMREPLICATE
;
; AUTHOR:
; Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;
; PURPOSE:
; Replicates an array or scalar into a larger array, as REPLICATE does.
;
; CALLING SEQUENCE:
; ARRAY = CMREPLICATE(VALUE, DIMS)
;
; DESCRIPTION:
;
; The CMREPLICATE function constructs an array, which is filled with
; the specified VALUE template. CMREPLICATE is very similar to the
; built-in IDL function REPLICATE. However there are two
; differences:
;
; * the VALUE can be either scalar or an ARRAY.
;
; * the dimensions are specified as a single vector rather than
; individual function arguments.
;
; For example, if VALUE is a 2x2 array, and DIMS is [3,4], then the
; resulting array will be 2x2x3x4.
;
; INPUTS:
;
; VALUE - a scalar or array template of any type, to be replicated.
; NOTE: These two calls do not produce the same result:
; ARRAY = CMREPLICATE( 1, DIMS)
; ARRAY = CMREPLICATE([1], DIMS)
; In the first case the output dimensions will be DIMS and
; in the second case the output dimensions will be 1xDIMS
; (except for structures). That is, a vector of length 1 is
; considered to be different from a scalar.
;
; DIMS - Dimensions of output array (which are combined with the
; dimensions of the input VALUE template). If DIMS is not
; specified then VALUE is returned unchanged.
;
; RETURNS:
; The resulting replicated array.
;
; EXAMPLE:
; x = [0,1,2]
; help, cmreplicate(x, [2,2])
; <Expression> INT = Array[3, 2, 2]
; Explanation: The 3-vector x is replicated 2x2 times.
;
; x = 5L
; help, cmreplicate(x, [2,2])
; <Expression> LONG = Array[2, 2]
; Explanation: The scalar x is replicated 2x2 times.
;
; SEE ALSO:
;
; REPLICATE
;
; MODIFICATION HISTORY:
; Written, CM, 11 Feb 2000
; Fixed case when ARRAY is array of structs, CM, 23 Feb 2000
; Apparently IDL 5.3 can't return from execute(). Fixed, CM, 24 Feb
; 2000
; Corrected small typos in documentation, CM, 22 Jun 2000
; Removed EXECUTE() call by using feature of REBIN() new in IDL 5.6,
; (thanks to Dick Jackson) CM, 24 Apr 2009
; Remove some legacy code no longer needed after above change
; (RETVAL variable no longer defined; thanks to A. van Engelen),
; CM, 08 Jul 2009
;
;-
; Copyright (C) 2000, 2009, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
function cmreplicate, array, dims
if n_params() EQ 0 then begin
message, 'RARRAY = CMREPLICATE(ARRAY, DIMS)', /info
return, 0L
endif
on_error, 2
if n_elements(dims) EQ 0 then return, array
if n_elements(array) EQ 0 then $
message, 'ERROR: ARRAY must have at least one element'
;; Construct new dimensions, being careful about scalars
sz = size(array)
type = sz(sz(0)+1)
if sz(0) EQ 0 then return, make_array(value=array, dimension=dims)
onedims = [sz(1:sz(0)), dims*0+1] ;; For REFORM, to extend # of dims.
newdims = [sz(1:sz(0)), dims ] ;; For REBIN, to enlarge # of dims.
nnewdims = n_elements(newdims)
if nnewdims GT 8 then $
message, 'ERROR: resulting array would have too many dimensions.'
if type NE 7 AND type NE 8 AND type NE 10 AND type NE 11 then begin
;; Handle numeric types
;; Passing REBIN(...,ARRAY) is a feature introduced in
;; IDL 5.6, and will be incompatible with previous versions.
array1 = rebin(reform([array], onedims), newdims)
return, array1
endif else begin
;; Handle strings, structures, pointers, and objects separately
;; Handle structures, which are never scalars
if type EQ 8 AND sz(0) EQ 1 AND n_elements(array) EQ 1 then $
return, reform(make_array(value=array, dimension=dims), dims, /over)
nold = n_elements(array)
nadd = 1L
for i = 0L, n_elements(dims)-1 do nadd = nadd * round(dims(i))
if type EQ 8 then $
array1 = make_array(value=array(0), dimension=[nold,nadd]) $
else $
array1 = make_array(type=type, dimension=[nold,nadd], /nozero)
array1 = reform(array1, [nold,nadd], /overwrite)
array2 = reform([array], n_elements(array))
;; Efficient copying, done by powers of two
array1(0,0) = temporary(array2)
stride = 1L ;; stride increase by a factor of two in each iteration
i = 1L & nleft = nadd - 1
while nleft GT stride do begin
array1(0,i) = array1(*,0:stride-1) ;; Note sneaky IDL optimization
i = i + stride & nleft = nleft - stride
stride = stride * 2
endwhile
if nleft GT 0 then array1(0,i) = array1(*,0:nleft-1)
return, reform(array1, newdims, /overwrite)
endelse
end