C +++ C C Source: src/lib/cubspl.F C C ---------------------------------------------- C SHADOW C Center for X-ray Lithography C University of Wisconsin-Madison C 3731 Schneider Dr., Stoughton, WI, 53589 C ---------------------------------------------- C C Log: cubspl.F C Revision 1.6 1991/07/06 19:56:48 khan C Grenoble and after. Minor changes C C Revision 1.5 91/04/05 13:54:17 cwelnak C changed quotes on #include C C Revision 1.4 91/03/22 12:08:43 cwelnak C SUN version -- INC to #inc C C Revision 1.3 90/11/13 14:04:48 khan C Cleanup and SAVE statements C C Revision 1.2 90/07/20 22:05:04 khan C put #if unix ... to make it work also on vms C C Revision 1.1 90/07/10 14:56:00 khan C Initial revision C C C --- #if defined(unix) || HAVE_F77_CPP # include #elif defined(vms) INCLUDE 'SHADOW$INC:HEADER.TXT/LIST' #endif C+++ C PROGRAM CUBSPL C C PURPOSE GENERATE THE COEFFICIENTS G(2,I), G(3,I), G(4,I) AND C G(5,I) FOR THE POLYNOMIAL Y=G(2,I)+X(I)*(G(3,I)+X(I) C *(G(4,I)+X(I)*G(5,I))) WHICH IS VALID BETWEEN THE C INTERVAL X(I) AND X(I+1). C C INPUT N =# OF PAIR OF RAW DATA POINTS (X(I),Y(I)). C Y(N) =THE ARRAY WHICH CONTAINS THE DATA OF Y(1) TO Y(N). C G(5,N)=THE FIRST ROW G(1,I) SHOULD CONTAIN THE DATA C OF X(1) TO X(N). THE OTHER FOUR ROWS WILL BE C FILLED UP WITH THE POLYNOMIAL COEFFICIENTS BY THIS C PROGRAM. C IER=1(0), check (no check) for steep slope at the two ends. C C OUTPUT THE LAST FOUR ROWS OF G. C IER=1 FOR ERROR, =0 OTHERWISE. C C--- SUBROUTINE CUBSPL(G, Y, N, IER) IMPLICIT REAL*8 (A-H,O-Z) #if defined(unix) || HAVE_F77_CPP # include #elif defined(vms) INCLUDE 'SHADOW$INC:DIM.PAR/LIST' INTEGER MPURGE(2) #endif REAL*8 G(5,N), Y(N), E1(N_DIM), E2(N_DIM), R INTEGER I, IER, N C IF (N.GT.N_DIM) CALL LEAVE $ ('CUBSPL','Splines cannot handle more than N_DIM points.' $ ,-1) IF (N.LT.4) CALL LEAVE $ ('CUBSPL','At least 4 data points are needed for splines.',-1) DO 11 J = 1, N-1 IF (G(1,J).GT.G(1,J+1)) THEN DO 21 I = 1, N E1(I) = G(1,I) E2(I) = Y(I) 21 CONTINUE CALL SORT_SPL (E1,E2,N) DO 31 I = 1, N G(1,I) = E1(I) Y(I) = E2(I) 31 CONTINUE GO TO 101 END IF 11 CONTINUE 101 SMIN = 1.0D+30 DO 41 I = 1, N-1 E1(I) = G(1,I+1) - G(1,I) C check for zero here ... IF (E1(I).EQ.0.0) THEN E1(I) = 1.0E-12 ENDIF C end changes. may need at a later date for special wigglers. clw C 7/22/93 E2(I) = (Y(I+1) - Y(I)) / E1(I) SMIN = MIN(SMIN,ABS(E2(I))) 41 CONTINUE IF (IER.EQ.0) THEN ISTART = 1 IEND = N ELSE C C Check if the slopes at the two ends (1 -> ISTART, IEND -> N) are too steep. C I = 1 51 IF (ABS(E2(I)).GT.SMIN*1) THEN I = I + 1 GOTO 51 END IF ISTART = I I = N - 1 61 IF (ABS(E2(I)).GT.SMIN*1) THEN I = I - 1 GOTO 61 END IF IEND = I IF ((IEND-ISTART+1).LT.4) THEN ISTART = 1 IEND =N END IF END IF C C Start computing the spline coefficients. C G(3,ISTART) = E1(ISTART+1) G(4,ISTART) = E1(ISTART) + E1(ISTART+1) G(5,ISTART) = ((E1(ISTART)+2.0D0*G(4,ISTART))*E1(ISTART+1)* $ E2(ISTART)+(E1(ISTART)**2.0D0) $ *E2(ISTART+1))/G(4,ISTART) DO 71 I = ISTART+1, IEND-1 G(2,I) = E1(I) G(3,I) = 2.0D0*(E1(I-1) + E1(I)) G(4,I) = E1(I-1) G(5,I) = 3.0D0*(E1(I)*E2(I-1) + E1(I-1)*E2(I)) 71 CONTINUE G(2,IEND) = E1(IEND-1) + E1(IEND-2) G(3,IEND) = E1(IEND-2) G(5,IEND) = ((E1(IEND-1)+2.0D0*G(2,IEND))*E2(IEND-1)*E1(IEND-2)+ $ (E1(IEND-1) $ **2.0D0)*E2(IEND-2))/G(2,IEND) DO 81 I = ISTART, IEND-1 R = G(2,I+1) / G(3,I) G(3,I+1) = G(3,I+1) - R*G(4,I) G(5,I+1) = G(5,I+1) - R*G(5,I) 81 CONTINUE G(3,IEND) = G(5,IEND) / G(3,IEND) DO 91 I = IEND-1, ISTART, -1 G(3,I) = (G(5,I) - G(4,I)*G(3,I+1)) / G(3,I) 91 CONTINUE DO 141 I = ISTART, IEND-1 G(2,I) = Y(I) G(4,I) = (3.0D0*E2(I) - 2.0D0*G(3,I) - G(3,I+1)) / E1(I) G(5,I) = (-2.0D0*E2(I) + G(3,I) + G(3,I+1)) / (E1(I)**2.0D0) 141 CONTINUE C C Use broken line instead at the two ends if slope are too steep. C DO 111 I = 1, ISTART-1 G(2,I) = Y(I) G(3,I) = E2(I) G(4,I) = 0.0D0 G(5,I) = 0.0D0 111 CONTINUE DO 121 I = IEND, N-1 G(2,I) = Y(I) G(3,I) = E2(I) G(4,I) = 0.0D0 G(5,I) = 0.0D0 121 CONTINUE IER = 0 #if defined(vms) MPURGE(1) = %LOC(E1(1)) MPURGE(2) = %LOC(E1(N_DIM)) CALL SYS$PURGWS (MPURGE) MPURGE(1) = %LOC(E2(1)) MPURGE(2) = %LOC(E2(N_DIM)) CALL SYS$PURGWS (MPURGE) #endif RETURN END ******************************************************************************* C+++ C SUBROUTINE SORT_SPL (Real*8 version) C C PURPOSE To sort a pair of array XVEC and YVEC. The result is C in ascending order according to XVEC. C--- SUBROUTINE SORT_SPL (XVEC,YVEC,ICOUNT) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION XVEC(ICOUNT),YVEC(ICOUNT) DO 11 I = 1, ICOUNT IMIN = I AMIN = XVEC(I) DO 21 J = I, ICOUNT IF (XVEC(J).LT.AMIN) THEN AMIN = XVEC(J) IMIN = J END IF 21 CONTINUE XTEMP = XVEC(I) XVEC(I) = XVEC(IMIN) XVEC(IMIN) = XTEMP YTEMP = YVEC(I) YVEC(I) = YVEC(IMIN) YVEC(IMIN) = YTEMP 11 CONTINUE RETURN END