aboutsummaryrefslogtreecommitdiff
path: root/src/xutils.f
blob: 3ad8d76ade76164136156f4ccdea80dc0d15d3a3 (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



      SUBROUTINE SETEXP(S,DS1,SMAX,NN)
C........................................................
C     Sets geometrically stretched array S:
C
C       S(i+1) - S(i)  =  r * [S(i) - S(i-1)]
C
C       S     (output)  array to be set  
C       DS1   (input)   first S increment:  S(2) - S(1)
C       SMAX  (input)   final S value:      S(NN)
C       NN    (input)   number of points
C........................................................
      REAL S(NN)
C
      SIGMA = SMAX/DS1
      NEX = NN-1
      RNEX = FLOAT(NEX)
      RNI = 1.0/RNEX
C
C---- solve quadratic for initial geometric ratio guess
      AAA = RNEX*(RNEX-1.0)*(RNEX-2.0) / 6.0
      BBB = RNEX*(RNEX-1.0) / 2.0
      CCC = RNEX - SIGMA
C
      DISC = BBB**2 - 4.0*AAA*CCC
      DISC = MAX( 0.0 , DISC )
C
      IF(NEX.LE.1) THEN
       STOP 'SETEXP: Cannot fill array.  N too small.'
      ELSE IF(NEX.EQ.2) THEN
       RATIO = -CCC/BBB  +  1.0
      ELSE
       RATIO = (-BBB + SQRT(DISC))/(2.0*AAA)  +  1.0
      ENDIF
C
      IF(RATIO.EQ.1.0) GO TO 11
C
C---- Newton iteration for actual geometric ratio
      DO 1 ITER=1, 100
        SIGMAN = (RATIO**NEX - 1.0) / (RATIO - 1.0)
        RES = SIGMAN**RNI - SIGMA**RNI
        DRESDR = RNI*SIGMAN**RNI
     &         * (RNEX*RATIO**(NEX-1) - SIGMAN) / (RATIO**NEX - 1.0)
C
        DRATIO = -RES/DRESDR
        RATIO = RATIO + DRATIO
C
        IF(ABS(DRATIO) .LT. 1.0E-5) GO TO 11
C
    1 CONTINUE
      WRITE(*,*) 'SETEXP: Convergence failed.  Continuing anyway ...'
C
C---- set up stretched array using converged geometric ratio
   11 S(1) = 0.0
      DS = DS1
      DO 2 N=2, NN
        S(N) = S(N-1) + DS
        DS = DS*RATIO
    2 CONTINUE
C
      RETURN
      END



      FUNCTION ATANC(Y,X,THOLD)
      IMPLICIT REAL (A-H,M,O-Z)
C---------------------------------------------------------------
C     ATAN2 function with branch cut checking.
C
C     Increments position angle of point X,Y from some previous
C     value THOLD due to a change in position, ensuring that the
C     position change does not cross the ATAN2 branch cut
C     (which is in the -x direction).  For example:
C
C       ATANC( -1.0 , -1.0 , 0.75*pi )  returns  1.25*pi , whereas
C       ATAN2( -1.0 , -1.0 )            returns  -.75*pi .
C
C     Typically, ATANC is used to fill an array of angles:
C
C        THETA(1) = ATAN2( Y(1) , X(1) )
C        DO i=2, N
C          THETA(i) = ATANC( Y(i) , X(i) , THETA(i-1) )
C        END DO
C
C     This will prevent the angle array THETA(i) from jumping by 
C     +/- 2 pi when the path X(i),Y(i) crosses the negative x axis.
C
C     Input:
C       X,Y     point position coordinates
C       THOLD   position angle of nearby point
C
C     Output:
C       ATANC   position angle of X,Y
C---------------------------------------------------------------
      DATA  PI /3.1415926535897932384/
      DATA TPI /6.2831853071795864769/
C
C---- set new position angle, ignoring branch cut in ATAN2 function for now
      THNEW = ATAN2( Y , X )
      DTHET = THNEW - THOLD
C
C---- angle change cannot exceed +/- pi, so get rid of any multiples of 2 pi 
      DTCORR = DTHET - TPI*INT( (DTHET + SIGN(PI,DTHET))/TPI )
C
C---- set correct new angle
      ATANC = THOLD + DTCORR
C
      RETURN
      END ! ATANC