aboutsummaryrefslogtreecommitdiff
path: root/src/sort.f
blob: a60ff1fac7017abf1a57e7217458edf53ae252f8 (plain)
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
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255

      SUBROUTINE HSORT(N,A,INDX)
      DIMENSION A(*)
      DIMENSION INDX(*)
C--------------------------------------
C     Heapsort algorithm.
C     Returns INDX(.) such that
C
C       A(INDX(i)) < A(INDX(i+1))
C
C     Stolen from Numerical Recipes.
C--------------------------------------
C
      DO I = 1, N
        INDX(I) = I
      ENDDO
C
      IF(N.LE.1) RETURN
C
      L = N/2 + 1
      IR = N
C
 10   CONTINUE
      IF(L.GT.1) THEN
        L = L-1
        INDXT = INDX(L)
        Q = A(INDXT)
      ELSE
        INDXT = INDX(IR)
        Q = A(INDXT)
        INDX(IR) = INDX(1)
C
        IR = IR - 1
        IF(IR.EQ.1) THEN
          INDX(1) = INDXT
          RETURN
        ENDIF
      ENDIF
C
      I = L
      J = L+L
C
 20   IF(J.LE.IR) THEN
        IF(J.LT.IR) THEN
          IF(A(INDX(J)) .LT. A(INDX(J+1))) J = J+1
        ENDIF
        IF(Q .LT. A(INDX(J))) THEN
          INDX(I) = INDX(J)
C
          I = J
          J = J+J
        ELSE
          J = IR+1
        ENDIF
        GO TO 20
      ENDIF
C
      INDX(I) = INDXT
      GO TO 10
      END

      SUBROUTINE ASORT(N,A,INDX,ATMP)
      DIMENSION A(*), ATMP(*)
      DIMENSION INDX(*)
C-----------------------------------------------
C     Applies sorted index array to reorder A.
C-----------------------------------------------
      DO I = 1, N
        ATMP(I) = A(I)
      ENDDO
C
      DO I = 1, N
        ISORT = INDX(I)
        A(I) = ATMP(ISORT)
      ENDDO
C
      RETURN
      END

      SUBROUTINE REMD(N,A,INDX,TOL,NNEW)
      DIMENSION A(*)
      DIMENSION INDX(*)
C----------------------------------------------------
C     Sets index array, such that 
C     duplicate A values are left out
C----------------------------------------------------
      K = 1
      INDX(K) = 1
C
      DO I = 2, N
        IF(ABS(A(I)-A(I-1)) .GT. TOL) THEN
          K = K + 1
          INDX(K) = I
        ENDIF
      ENDDO
C
      NNEW = K
C
      RETURN
      END ! REMD


      SUBROUTINE SORTDUP(KK,S,W)
C--- Sort arrays in S with no removal of duplicates
      DIMENSION S(KK), W(KK)
      LOGICAL DONE
C
C---- sort arrays
      DO 10 IPASS=1, 1234
        DONE = .TRUE.
        DO 101 N=1, KK-1
          NP = N+1
          IF(S(NP).GE.S(N)) GO TO 101
           TEMP = S(NP)
           S(NP) = S(N)
           S(N) = TEMP
           TEMP = W(NP)
           W(NP) = W(N)
           W(N) = TEMP
           DONE = .FALSE.
  101   CONTINUE
        IF(DONE) GO TO 11
   10 CONTINUE
      WRITE(*,*) 'Sort failed'
C
   11 CONTINUE
      RETURN
      END


      SUBROUTINE FIXDUP(KK,S,W)
C--- Check arrays in S by removing leading and ending duplicates
C    eliminate extra duplicates (more than one duplicate point) elsewhere
      DIMENSION S(KK), W(KK)
      LOGICAL DONE
C
C---- Check first elements for dups
      IF(S(2).EQ.S(1)) THEN
        DO N=1, KK-1
          S(N) = S(N+1)
          W(N) = W(N+1)
        END DO
        KK = KK - 1
      ENDIF
C
C---- Check last elements for dups
      IF(S(KK).EQ.S(KK-1)) THEN
        S(KK-1) = S(KK)
        W(KK-1) = W(KK)
        KK = KK - 1
      ENDIF
C
C--- Eliminate more than 2 succeeding identical elements 
 10   DO N=1, KK-2
        IF(S(N).EQ.S(N+1) .AND. S(N).EQ.S(N+2)) THEN
          DO I = N, KK-1
           S(I) = S(I+1)
           W(I) = W(I+1)
          END DO
          KK = KK - 1
          GO TO 10
        ENDIF
      END DO
C
      RETURN
      END


      SUBROUTINE SORT(KK,S,W)
      DIMENSION S(KK), W(KK)
      LOGICAL DONE
C
C---- sort arrays
      DO 10 IPASS=1, 1234
        DONE = .TRUE.
        DO 101 N=1, KK-1
          NP = N+1
          IF(S(NP).GE.S(N)) GO TO 101
           TEMP = S(NP)
           S(NP) = S(N)
           S(N) = TEMP
           TEMP = W(NP)
           W(NP) = W(N)
           W(N) = TEMP
           DONE = .FALSE.
  101   CONTINUE
        IF(DONE) GO TO 11
   10 CONTINUE
      WRITE(*,*) 'Sort failed'
C
C---- search for duplicate pairs and eliminate each one
   11 KKS = KK
      DO 20 K=1, KKS
        IF(K.GE.KK) RETURN
        IF(S(K).NE.S(K+1)) GO TO 20
C------- eliminate pair
         KK = KK-2
         DO 201 KT=K, KK
           S(KT) = S(KT+2)
           W(KT) = W(KT+2)
  201    CONTINUE
   20 CONTINUE
C
      RETURN
      END


 
      SUBROUTINE SORTOL(TOL,KK,S,W)
      DIMENSION S(KK), W(KK)
      LOGICAL DONE
C
C---- sort arrays
      DO IPASS=1, 1234
        DONE = .TRUE.
        DO N=1, KK-1
          NP = N+1
          IF(S(NP).LT.S(N)) THEN
           TEMP = S(NP)
           S(NP) = S(N)
           S(N) = TEMP
           TEMP = W(NP)
           W(NP) = W(N)
           W(N) = TEMP
           DONE = .FALSE.
          ENDIF
        END DO
        IF(DONE) GO TO 10
      END DO
      WRITE(*,*) 'Sort failed'
C
C---- search for near-duplicate pairs and eliminate extra points
C---- Modified 4/24/01 HHY to check list until ALL duplicates removed
C     This cures a bug for sharp LE foils where there were 3 LE points in
C     camber, thickness lists from GETCAM.
C
 10   KKS = KK
      DONE = .TRUE.
      DO 20 K=1, KKS
        IF(K.GE.KK) GO TO 20
        DSQ = (S(K)-S(K+1))**2 + (W(K)-W(K+1))**2
        IF(DSQ.GE.TOL*TOL) GO TO 20
C------- eliminate extra point pairs
ccc         write(*,*) 'extra on point ',k,kks
         KK = KK-1
         DO KT=K+1, KK
           S(KT) = S(KT+1)
           W(KT) = W(KT+1)
         END DO
         DONE = .FALSE.
   20 CONTINUE
      IF(.NOT.DONE) GO TO 10
C
      RETURN
      END