C @(#)necwlcopt.for 17.1.1.1 (ES0-DMD) 01/25/02 17:51:32 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.COPYRIGHT (C) 1992 European Southern Observatory C.IDENT wlcoptim.for C.AUTHOR Pascal Ballester, ESO - Garching C.KEYWORDS Spectroscopy, Echelle, C.PURPOSE Computes optimal parameters for the C wavelength calibration C.VERSION 1.0 Creation 1992-MAY-26 C------------------------------------------------------- C C C Program : Wavelength Calibration Least Squares Solutions C C This program has been partly generated using Macsyma to C solve the following least-square problem: C C Xr(i) = cos(a) * X(i) + bin * sin(a) * Y(i) (1) C C takes into account a binning factor of the ccd and a rotation C of the detector. C C F(Xr(i)) = B * Xr(i) + A = (m+p(i))*lambda(i) (2) C C is the echelle relation applied in the rotated space and approximated C by a linear relation. m is the absolute order number of a line C and p(i) is relative order number of the other lines. C C Macsyma has been used to find different optima, like: C - The absolute order number of the first identification : m C - The binning factor : bin C - In method ANGLE, the angle of rotation : a C - In method PAIR, the wavelength of the two pairs : lambda(1), lambda(2) C C The different least-square solutions are used to compute the optimal C value of a parameter, assuming that all other parameters are correctly C provided. This is a cross-check of the entered values. By computing C the new rms after replacement of a parameter by its optimal value, C it is possible to provide a diagnosis of what action would benefit C the most to solve the problem. C C INPUT: C IN_A/C/1/60 Name of line table C IN_B/C/1/60 Method (PAIR or ANGLE) C INPUTR/R/1/2 Angle of rotation, Binning ratio. C C OUTPUT: C OUTPUTI/I/1/1 Order number of the first line C OUTPUTR/R/1/2 Optimal rotation angle, Optimal binning factor C IMPLICIT NONE INTEGER MAXDIM PARAMETER (MAXDIM=1500) DOUBLE PRECISION IX(MAXDIM), IY(MAXDIM) DOUBLE PRECISION IP(MAXDIM), IL(MAXDIM) DOUBLE PRECISION SX, SY, SX2, SY2, SXY C DOUBLE PRECISION SML2K, SXMLK, SYMLK, SMLK DOUBLE PRECISION SP2L2K, BUFFER DOUBLE PRECISION SXLK, SXPLK, SYLK, SYPLK, SPLK, SLK DOUBLE PRECISION SL2K, SPL2K, PK, LK, XK, YK, RMSLK DOUBLE PRECISION SML2, SP2L2, SXML, SYML, SML DOUBLE PRECISION SXL, SXPL, SYL, SYPL, SPL, SL DOUBLE PRECISION SL2, SPL2, SBIN INTEGER M, N, K, M2OPT, LOOP, ERROR, P DOUBLE PRECISION BIN, PI, MDP, MINRMS(MAXDIM) DOUBLE PRECISION A, B1 DOUBLE PRECISION TGA1, TGA2, AOPT REAL OUTANG, RBIN DOUBLE PRECISION L1,L2,X1,X2,XP1,XP2,L1OPT,L2OPT DOUBLE PRECISION SXR, SXR2, SXRML, DET, COEFA, COEFB, RMS DOUBLE PRECISION L1OLD, L2OLD, RMSL1, RMSL2 CHARACTER MODE*5, METHOD*10, TEXT*80 INTEGER MOLD, ROW DOUBLE PRECISION RMSM, MRMS, ARMS, BRMS, AOLD, BOLD DOUBLE PRECISION LKOPT(MAXDIM), LKO CHARACTER TABLE*60 INTEGER IAV, KUN, KNUL, ISTAT, ICOLX, ICOLW INTEGER ICOLO, ICOLN, TID, NCOL1, NROW1, NSC1 INTEGER NACOL, NAROW, MADRID LOGICAL SELECT INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON/VMR/MADRID(1) INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C Initialisations C C ... Init input table C CALL STSPRO('WLCLS') CALL STKRDC('IN_A',1,1,60,IAV,TABLE,KUN,KNUL,ISTAT) CALL TBTOPN(TABLE,F_U_MODE,TID,ISTAT) C CALL TBIGET(TID,NCOL1,NROW1,NSC1,NACOL,NAROW,ISTAT) CALL TBCSER(TID,':X ' ,ICOLX,ISTAT) CALL TBCSER(TID,':IDENT' ,ICOLW,ISTAT) CALL TBCSER(TID,':ORDER' ,ICOLO,ISTAT) CALL TBCSER(TID,':YNEW' ,ICOLN,ISTAT) CALL STKRDC('IN_B',1,1,10,IAV,METHOD,KUN,KNUL,ISTAT) CALL STKRDR('INPUTR',1,1,IAV,OUTANG,KUN,KNUL,ISTAT) ! Rotation angle CALL STKRDR('INPUTR',2,1,IAV,RBIN,KUN,KNUL,ISTAT) ! Binning factor PI = 3.14159 A = OUTANG*PI/180. BIN = RBIN ERROR = 0 IF (METHOD(1:1) .EQ. 'p') WRITE(METHOD,111) 'PAIR' IF (METHOD(1:1) .EQ. 'a') WRITE(METHOD,111) 'ANGLE' 111 FORMAT(A5) N = 0 DO 100 ROW = NROW1, 1, -1 CALL TBSGET(TID, ROW, SELECT, ISTAT) IF (SELECT) THEN N = N + 1 CALL TBERDD(TID, ROW, ICOLO, IP(N), KNUL, ISTAT) CALL TBERDD(TID, ROW, ICOLX, IX(N), KNUL, ISTAT) CALL TBERDD(TID, ROW, ICOLN, IY(N), KNUL, ISTAT) CALL TBERDD(TID, ROW, ICOLW, IL(N), KNUL, ISTAT) ENDIF 100 CONTINUE M = IP(1) DO 110 ROW = 1, N IP(ROW) = IP(ROW) - M 110 CONTINUE C Check maximum size of arrays IF (MAXDIM .LE. (N+10)) THEN CALL STTPUT( + 'Buffer size MAXDIM too small in module wlcoptim',ISTAT) GOTO 1900 ENDIF C Check organisation of data if method PAIR IF (METHOD(1:1) .EQ. 'P') THEN IF (IL(1) .NE. IL(2)) ERROR = 1 IF (IL(3) .NE. IL(4)) ERROR = 1 IF (IP(2) .NE. (IP(1)+1)) ERROR = 1 IF (IP(4) .NE. (IP(3)+1)) ERROR = 1 IF (N .NE. 4) ERROR = 1 IF (ERROR .EQ. 1) THEN IF (IL(1).EQ.IL(3).AND.IL(2).EQ.IL(4).AND.IP(3).EQ. C (IP(1)+1).AND.IP(4).EQ.(IP(2)+1).AND.N.EQ.4) THEN ERROR = 0 BUFFER = IL(2) IL(2) = IL(3) IL(3) = BUFFER BUFFER = IX(2) IX(2) = IX(3) IX(3) = BUFFER BUFFER = IY(2) IY(2) = IY(3) IY(3) = BUFFER BUFFER = IP(2) IP(2) = IP(3) IP(3) = NINT(BUFFER) ELSE CALL STTPUT( + 'Error: Identifications are wrongly organized',ISTAT) CALL STTPUT( + 'Remember : First line: left-hand side, then right-side', + ISTAT) CALL STTPUT( + ' Second line: left-hand side, then right-side', + ISTAT) GOTO 1900 ENDIF ENDIF ENDIF C Check that identifications are in different orders for method ANGLE IF (METHOD(1:1) .EQ. 'A') THEN SPL = 0.D0 DO 200 ROW = 1,N SPL = SPL + IP(ROW) 200 CONTINUE C234567890123456789012345678901234567890123456789012345678901234567890123456 IF (SPL .LT. 0.5) THEN CALL STTPUT( + 'Error: Identifications cannot be all in the same order',ISTAT) GOTO 1900 ENDIF ENDIF C Computes sums SX = 0.D0 SY = 0.D0 SX2 = 0.D0 SXY = 0.D0 SY2 = 0.D0 SYL = 0.D0 SYPL = 0.D0 SPL = 0.D0 SL = 0.D0 SL2 = 0.D0 SPL2 = 0.D0 SP2L2 = 0.D0 SXL = 0.D0 SXPL = 0.D0 DO 10 LOOP = 1 , N SX = SX + IX(LOOP) SY = SY + IY(LOOP) SX2 = SX2 + IX(LOOP)*IX(LOOP) SY2 = SY2 + IY(LOOP)*IY(LOOP) SXY = SXY + IX(LOOP)*IY(LOOP) SYL = SYL + IY(LOOP)*IL(LOOP) SYPL = SYPL + IY(LOOP)*IP(LOOP)*IL(LOOP) SPL = SPL + IP(LOOP)*IL(LOOP) SL = SL + IL(LOOP) SL2 = SL2 + IL(LOOP)*IL(LOOP) SPL2 = SPL2 + IP(LOOP)*IL(LOOP)**2 SP2L2 = SP2L2 + IP(LOOP)**2*IL(LOOP)**2 SXL = SXL + IX(LOOP)*IL(LOOP) SXPL = SXPL + IX(LOOP)*IP(LOOP)*IL(LOOP) 10 CONTINUE SXML = M*SXL + SXPL SYML = M*SYL + SYPL SML = M*SL + SPL SML2 = M*M*SL2 + 2*M*SPL2 + SP2L2 C Computes a linear fit and display results SXR = COS(A)*SX + BIN*SIN(A)*SY SXR2 = SX2*COS(A)**2+2*SXY*COS(A)*BIN*SIN(A)+ 1 SY2*SIN(A)**2*BIN**2 SXRML = SXML*COS(A)+SYML*SIN(A)*BIN DET = N*SXR2-SXR*SXR COEFA = (SML*SXR2-SXR*SXRML)/DET COEFB = (N*SXRML-SXR*SML)/DET RMS = SML2-2*COEFA*COEFB*SXR-COEFB*COEFB*SXR2-N*COEFA*COEFA RMS = SQRT(RMS/N) C Display header WRITE(TEXT,510) CALL STTPUT(TEXT,ISTAT) WRITE(TEXT,570) CALL STTPUT(TEXT,ISTAT) WRITE(TEXT,510) CALL STTPUT(TEXT,ISTAT) WRITE(TEXT,520) CALL STTPUT(TEXT,ISTAT) WRITE(TEXT,500) COEFA,COEFB,RMS CALL STTPUT(TEXT,ISTAT) WRITE(TEXT,510) CALL STTPUT(TEXT,ISTAT) WRITE(TEXT,530) CALL STTPUT(TEXT,ISTAT) WRITE(TEXT,510) CALL STTPUT(TEXT,ISTAT) 500 FORMAT (1X,'* A = ',F12.5,' B= ',F12.5,' RMS = ',F12.3,7X,'*') 510 FORMAT (1X,61('*')) 520 FORMAT (1X,'*',10X,'Linear fit of : m * lambda = b * xr + a', C 10X,'*') 530 FORMAT (1X,'* Parameter * Value * Optimum *', C ' RMS *') 540 FORMAT (1X,'*',A16,' *',F12.5,' *',F13.5,' *',F11.2,' *') 550 FORMAT (1X,'*',' Ident ',I6,4X,'*',F12.5,' *',F13.5,' *', C F11.2,' *') 560 FORMAT (1X,'*',A16,' * ',I6,' * ',I7,' *',F11.2,' *') 570 FORMAT (1X,'*',13X,'Wavelength Calibration Diagnosis',14X,'*') C Least square solutions for the methods PAIR and ANGLE C Optimal Rotation Angle DET = ABS(N*SYML-SML*SY) IF (DET .GT. 0.01) THEN TGA1 = (N*(SX2*SYML-SXML*SXY)-SX**2*SYML+SML*(SX*SXY-SX2 1 *SY)+SX*SXML*SY)/(BIN*(N*(SXML*SY2-SXY*SYML)+SX*SY*SYML+SML*(SX 2 Y*SY-SX*SY2)-SXML*SY**2)) TGA2 = -(N*SXML-SML*SX)/(BIN*(N*SYML-SML*SY)) AOPT = ATAN(TGA1) ELSE AOPT = -999 ENDIF C Optimal Binning Factor IF (ABS(A) .GT. 0.02) THEN B1 = COS(A)*(N*(SX2*SYML-SXML*SXY)-SX**2*SYML+SML*(SX*SXY-SX2*SY 1 )+SX*SXML*SY)/(SIN(A)*(N*(SXML*SY2-SXY*SYML)+SX*SY*SYML+SML*(SX 2 Y*SY-SX*SY2)-SXML*SY**2)) ELSE B1 = -999 ENDIF C Optimal Absolute Order Number for the first identified line MDP = -((SIN(A)**2*BIN**2*N*SYL-SIN(A)**2*BIN**2*SL*SY+COS( 1 A)*SIN(A)*BIN*N*SXL-COS(A)*SIN(A)*BIN*SL*SX)*SYPL+(-SIN(A)**2*B 2 IN**2*SPL*SY+COS(A)*SIN(A)*BIN*N*SXPL-COS(A)*SIN(A)*BIN*SPL*SX) 3 *SYL+(SIN(A)**2*BIN**2*SL*SPL-SIN(A)**2*BIN**2*N*SPL2)*SY2+SIN( 4 A)**2*BIN**2*SPL2*SY**2+(-COS(A)*SIN(A)*BIN*SL*SXPL-COS(A)*SIN( 5 A)*BIN*SPL*SXL+2*COS(A)*SIN(A)*BIN*SPL2*SX)*SY+(2*COS(A)*SIN(A) 6 *BIN*SL*SPL-2*COS(A)*SIN(A)*BIN*N*SPL2)*SXY+(COS(A)**2*N*SXL-CO 7 S(A)**2*SL*SX)*SXPL-COS(A)**2*SPL*SX*SXL+(COS(A)**2*SL*SPL-COS( 8 A)**2*N*SPL2)*SX2+COS(A)**2*SPL2*SX**2)/(SIN(A)**2*BIN**2*N*SYL 9 **2+(-2*SIN(A)**2*BIN**2*SL*SY+2*COS(A)*SIN(A)*BIN*N*SXL-2*COS( : A)*SIN(A)*BIN*SL*SX)*SYL+(SIN(A)**2*BIN**2*SL**2-SIN(A)**2*BIN* ; *2*N*SL2)*SY2+SIN(A)**2*BIN**2*SL2*SY**2+(2*COS(A)*SIN(A)*BIN*S < L2*SX-2*COS(A)*SIN(A)*BIN*SL*SXL)*SY+(2*COS(A)*SIN(A)*BIN*SL**2 = -2*COS(A)*SIN(A)*BIN*N*SL2)*SXY+COS(A)**2*N*SXL**2-2*COS(A)**2* > SL*SX*SXL+(COS(A)**2*SL**2-COS(A)**2*N*SL2)*SX2+COS(A)**2*SL2*S ? X**2) M2OPT = NINT(MDP) C Compute new rms if a value is changed. MOLD = M ! Keep the old value in memory M = M2OPT ! Replace current value by optimum MODE = 'M' ! Set the mode GOTO 1000 ! Compute rms 1001 CONTINUE ! Return point M = MOLD ! Restore old value MRMS = RMSM ! Store rms MINRMS(1) = RMSM WRITE(TEXT,560) 'Order number',M,M2OPT,MRMS CALL STTPUT(TEXT,ISTAT) CALL STKWRI('OUTPUTI',M2OPT,1,1,KUN,ISTAT) AOLD = A A = AOPT MODE = 'A' GOTO 1000 1002 CONTINUE ! Return point A = AOLD ARMS = RMSM MINRMS(2) = RMSM WRITE(TEXT,540) 'Rotation angle',(A*180./PI),(AOPT*180/PI),ARMS CALL STTPUT(TEXT,ISTAT) OUTANG = AOPT*180./PI CALL STKWRR('OUTPUTR',OUTANG,1,1,KUN,ISTAT) BOLD = BIN BIN = B1 MODE = 'B' GOTO 1000 1003 CONTINUE ! Return point BIN = BOLD BRMS = RMSM MINRMS(3) = RMSM IF (B1 .GT. 0) THEN WRITE(TEXT,540) 'Binning factor',BIN,SQRT(B1),BRMS CALL STTPUT(TEXT,ISTAT) ELSE WRITE(TEXT,540) 'Binning factor',BIN,BIN,0. CALL STTPUT(TEXT,ISTAT) ENDIF C C Least square solutions for the method PAIR IF (METHOD(1:1) .EQ. 'P') THEN X1 = IX(1) XP1 = IX(2) X2 = IX(3) XP2 = IX(4) P = IP(3) - IP(1) L1 = IL(1) L2 = IL(3) L1OPT = (((2*L2*M+L2)*P+2*L2*M**2+L2*M)*XP2**2+((L2*P+2*L2*M+2*L2) 1 *XP1+((-4*L2*M-2*L2)*P-4*L2*M**2-4*L2*M-L2)*X2+(-L2*P-L2)*X1) 2 *XP2+(2*L2*M*P+2*L2*M**2+L2*M)*XP1**2+((L2*P-L2)*X2+ 3 ((-4*L2*M-2*L2)*P-4*L2*M**2-4*L2*M-L2)*X1)*XP1+((2*L2*M+L2)*P+ 4 2*L2*M**2+3*L2*M+L2)*X2**2+(-L2*P-2*L2*M)*X1*X2+((2*L2*M+2*L2) 5 *P+2*L2*M**2+3*L2*M+L2)*X1**2)/((2*M**2+2*M+2)*XP2**2+(2*M*XP1 6 +(-4*M**2-4*M-2)*X2+(-2*M-2)*X1)*XP2+2*M**2*XP1**2+(2*M*X2+(-4 7 *M**2-4*M)*X1)*XP1+(2*M**2+2*M+2)*X2**2+(-2*M-2)*X1*X2+(2*M**2 8 +4*M+2)*X1**2) L1OLD = L1 L1 = L1OPT RMSL1 = -4*((2*L2*P+2*(L2+L1)*M+L2+L1)*(XP2**2+XP1**2+X2**2+X1**2) 1 -(XP2+XP1+X2+X1)*(L2*(P+M+1)*XP2+L1*(M+1)*XP1+L2*(P+M)*X2+L1*M*X1 2 ))**2/(4*(XP2**2+XP1**2+X2**2+X1**2)-(XP2+XP1+X2+X1)**2)**2-2*( 3 XP2+XP1+X2+X1)*(4*(L2*(P+M+1)*XP2+L1*(M+1)*XP1+L2*(P+M)*X2+L1*M 4 *X1)-(2*L2*P+2*(L2+L1)*M+L2+L1)*(XP2+XP1+X2+X1))*((2*L2*P+2*(L2 5 +L1)*M+L2+L1)*(XP2**2+XP1**2+X2**2+X1**2)-(XP2+XP1+X2+X1)*(L2*( 6 P+M+1)*XP2+L1*(M+1)*XP1+L2*(P+M)*X2+L1*M*X1))/(4*(XP2**2+XP1**2 7 +X2**2+X1**2)-(XP2+XP1+X2+X1)**2)**2-(4*(L2*(P+M+1)*XP2+L1*(M+1 8 )*XP1+L2*(P+M)*X2+L1*M*X1)-(2*L2*P+2*(L2+L1)*M+L2+L1)*(XP2+XP1+ 9 X2+X1))**2*(XP2**2+XP1**2+X2**2+X1**2)/(4*(XP2**2+XP1**2+X2**2+ : X1**2)-(XP2+XP1+X2+X1)**2)**2+L2**2*(P+M+1)**2+L2**2*(P+M)**2+L ; 1**2*(M+1)**2+L1**2*M**2 L1 = L1OLD L2OPT = (((2*L1*M+L1)*P+2*L1*M**2+L1*M)*XP2**2+((L1*P+2*L1*M+2 1 *L1)*XP1+((-4*L1*M-2*L1)*P-4*L1*M**2-4*L1*M-L1)*X2+(-L1*P-L1)*X 2 1)*XP2+(2*L1*M*P+2*L1*M**2+L1*M)*XP1**2+((L1*P-L1)*X2+((-4*L1*M 3 -2*L1)*P-4*L1*M**2-4*L1*M-L1)*X1)*XP1+((2*L1*M+L1)*P+2*L1*M**2+ 4 3*L1*M+L1)*X2**2+(-L1*P-2*L1*M)*X1*X2+((2*L1*M+2*L1)*P+2*L1*M** 5 2+3*L1*M+L1)*X1**2)/((2*P**2+4*M*P+2*M**2)*XP2**2+((2*P+2*M)*XP 6 1+(-4*P**2+(-8*M-4)*P-4*M**2-4*M)*X2+(2*P+2*M)*X1)*XP2+(2*P**2+ 7 (4*M+2)*P+2*M**2+2*M+2)*XP1**2+((-2*P-2*M-2)*X2+(-4*P**2+(-8*M- 8 4)*P-4*M**2-4*M-2)*X1)*XP1+(2*P**2+(4*M+4)*P+2*M**2+4*M+2)*X2** 9 2+(-2*P-2*M-2)*X1*X2+(2*P**2+(4*M+2)*P+2*M**2+2*M+2)*X1**2) L2OLD = L2 L2 = L2OPT RMSL2 = -4*((2*L2*P+2*(L2+L1)*M+L2+L1)*(XP2**2+XP1**2+X2**2+X1**2) 1 -(XP2+XP1+X2+X1)*(L2*(P+M+1)*XP2+L1*(M+1)*XP1+L2*(P+M)*X2+L1*M*X1 2 ))**2/(4*(XP2**2+XP1**2+X2**2+X1**2)-(XP2+XP1+X2+X1)**2)**2-2*( 3 XP2+XP1+X2+X1)*(4*(L2*(P+M+1)*XP2+L1*(M+1)*XP1+L2*(P+M)*X2+L1*M 4 *X1)-(2*L2*P+2*(L2+L1)*M+L2+L1)*(XP2+XP1+X2+X1))*((2*L2*P+2*(L2 5 +L1)*M+L2+L1)*(XP2**2+XP1**2+X2**2+X1**2)-(XP2+XP1+X2+X1)*(L2*( 6 P+M+1)*XP2+L1*(M+1)*XP1+L2*(P+M)*X2+L1*M*X1))/(4*(XP2**2+XP1**2 7 +X2**2+X1**2)-(XP2+XP1+X2+X1)**2)**2-(4*(L2*(P+M+1)*XP2+L1*(M+1 8 )*XP1+L2*(P+M)*X2+L1*M*X1)-(2*L2*P+2*(L2+L1)*M+L2+L1)*(XP2+XP1+ 9 X2+X1))**2*(XP2**2+XP1**2+X2**2+X1**2)/(4*(XP2**2+XP1**2+X2**2+ : X1**2)-(XP2+XP1+X2+X1)**2)**2+L2**2*(P+M+1)**2+L2**2*(P+M)**2+L ; 1**2*(M+1)**2+L1**2*M**2 RMSL1 = SQRT(RMSL1/4.) RMSL2 = SQRT(RMSL2/4.) WRITE(TEXT,540) 'Ident. 1 and 2',L1OLD,L1OPT,RMSL1 CALL STTPUT(TEXT,ISTAT) WRITE(TEXT,540) 'Ident. 3 and 4',L2OLD,L2OPT,RMSL2 CALL STTPUT(TEXT,ISTAT) MINRMS(4) = RMSL1 MINRMS(5) = RMSL1 MINRMS(6) = RMSL2 MINRMS(7) = RMSL2 LKOPT(1) = L1OPT LKOPT(2) = L1OPT LKOPT(3) = L2OPT LKOPT(4) = L2OPT ENDIF C Optimal Wavelength for each of the computed wavelength IL(LOOP) IF (METHOD(1:1) .EQ. 'A') THEN DO 50 K = 1,N PK = IP(K) LK = IL(K) XK = IX(K) YK = IY(K) SYLK = 0.D0 SYPLK = 0.D0 SPLK = 0.D0 SLK = 0.D0 SL2K = 0.D0 SPL2K = 0.D0 SP2L2K = 0.D0 SXLK = 0.D0 SXPLK = 0.D0 DO 30 LOOP = 1 , N IF (LOOP.NE.K) THEN SYLK = SYLK + IY(LOOP)*IL(LOOP) SYPLK = SYPLK + IY(LOOP)*IP(LOOP)*IL(LOOP) SPLK = SPLK + IP(LOOP)*IL(LOOP) SLK = SLK + IL(LOOP) SL2K = SL2K + IL(LOOP)*IL(LOOP) SPL2K = SPL2K + IP(LOOP)*IL(LOOP)**2 SP2L2K = SP2L2K + IP(LOOP)**2*IL(LOOP)**2 SXLK = SXLK + IX(LOOP)*IL(LOOP) SXPLK = SXPLK + IX(LOOP)*IP(LOOP)*IL(LOOP) ENDIF 30 CONTINUE C The formula for LKOPT is splitted in two parts (too big for the poor C compiler !!) SBIN = ( < (SIN(A)**2*BIN**2*N*PK+SIN(A)**2*BIN**2*M*N)*YK**2+((2*COS(A)*S = IN(A)*BIN*N*PK+2*COS(A)*SIN(A)*BIN*M*N)*XK+(-2*SIN(A)**2*BIN**2 > *PK-2*SIN(A)**2*BIN**2*M)*SY+(-2*COS(A)*SIN(A)*BIN*PK-2*COS(A)* ? SIN(A)*BIN*M)*SX)*YK+(COS(A)**2*N*PK+COS(A)**2*M*N)*XK**2+((-2* @ COS(A)*SIN(A)*BIN*PK-2*COS(A)*SIN(A)*BIN*M)*SY+(-2*COS(A)**2*PK 1 -2*COS(A)**2*M)*SX)*XK+((SIN(A)**2*BIN**2-SIN(A)**2*BIN**2*N)*P 2 K-SIN(A)**2*BIN**2*M*N+SIN(A)**2*BIN**2*M)*SY2+(SIN(A)**2*BIN** 3 2*PK+SIN(A)**2*BIN**2*M)*SY**2+(2*COS(A)*SIN(A)*BIN*PK+2*COS(A) 4 *SIN(A)*BIN*M)*SX*SY+((2*COS(A)*SIN(A)*BIN-2*COS(A)*SIN(A)*BIN* 5 N)*PK-2*COS(A)*SIN(A)*BIN*M*N+2*COS(A)*SIN(A)*BIN*M)*SXY+((COS( 6 A)**2-COS(A)**2*N)*PK-COS(A)**2*M*N+COS(A)**2*M)*SX2+(COS(A)**2 7 *PK+COS(A)**2*M)*SX**2) LKOPT(K) = -((SIN(A)**2*BIN**2*N*SYPLK+SIN(A)**2*BIN**2*M*N*SY 1 LK+(-SIN(A)**2*BIN**2*SPLK-SIN(A)**2*BIN**2*M*SLK)*SY+COS(A)*SI 2 N(A)*BIN*N*SXPLK+COS(A)*SIN(A)*BIN*M*N*SXLK+(-COS(A)*SIN(A)*BIN 3 *SPLK-COS(A)*SIN(A)*BIN*M*SLK)*SX)*YK+(COS(A)*SIN(A)*BIN*N*SYPL 4 K+COS(A)*SIN(A)*BIN*M*N*SYLK+(-COS(A)*SIN(A)*BIN*SPLK-COS(A)*SI 5 N(A)*BIN*M*SLK)*SY+COS(A)**2*N*SXPLK+COS(A)**2*M*N*SXLK+(-COS(A 6 )**2*SPLK-COS(A)**2*M*SLK)*SX)*XK+(-SIN(A)**2*BIN**2*SY-COS(A)* 7 SIN(A)*BIN*SX)*SYPLK+(-SIN(A)**2*BIN**2*M*SY-COS(A)*SIN(A)*BIN* 8 M*SX)*SYLK+(SIN(A)**2*BIN**2*SPLK+SIN(A)**2*BIN**2*M*SLK)*SY2+( 9 -COS(A)*SIN(A)*BIN*SXPLK-COS(A)*SIN(A)*BIN*M*SXLK)*SY+(2*COS(A) : *SIN(A)*BIN*SPLK+2*COS(A)*SIN(A)*BIN*M*SLK)*SXY-COS(A)**2*SX*SX ; PLK-COS(A)**2*M*SX*SXLK+(COS(A)**2*SPLK+COS(A)**2*M*SLK)*SX2)/ + SBIN LKO = LKOPT(K) RMSLK = -N*((SPLK+M*SLK+LKO*(PK+M))*(SIN(A)**2*BIN**2*SY2+2*COS(A) 1 *SIN(A)*BIN*SXY+COS(A)**2*SX2)-(SIN(A)*BIN*SY+COS(A)*SX)*(SIN(A 2 )*BIN*(M*(LKO*YK+SYLK)+LKO*PK*YK+SYPLK)+COS(A)*(M*(LKO*XK+SXLK) C +LKO 3 *PK*XK+SXPLK)))**2/(N*(SIN(A)**2*BIN**2*SY2+2*COS(A)*SIN(A)*BIN 4 *SXY+COS(A)**2*SX2)-(SIN(A)*BIN*SY+COS(A)*SX)**2)**2-2*(SIN(A)* 5 BIN*SY+COS(A)*SX)*(N*(SIN(A)*BIN*(M*(LKO*YK+SYLK)+LKO*PK*YK+SYPLK 6)+COS(A)*(M*(LKO*XK+SXLK)+LKO*PK*XK+SXPLK))-(SPLK+M*SLK+LKO*(PK+M) 7 )*(SIN(A)*BIN*SY+COS(A)*SX))*((SPLK+M*SLK+LKO*(PK+M))*(SIN(A)**2 8 *BIN**2*SY2+2*COS(A)*SIN(A)*BIN*SXY+COS(A)**2*SX2)-(SIN(A)*BIN* 9 SY+COS(A)*SX)*(SIN(A)*BIN*(M*(LKO*YK+SYLK)+LKO*PK*YK+SYPLK)+COS(A : )*(M*(LKO*XK+SXLK)+LKO*PK*XK+SXPLK)))/(N*(SIN(A)**2*BIN**2*SY2+2* ; COS(A)*SIN(A)*BIN*SXY+COS(A)**2*SX2)-(SIN(A)*BIN*SY+COS(A)*SX)* < *2)**2-(SIN(A)**2*BIN**2*SY2+2*COS(A)*SIN(A)*BIN*SXY+COS(A)**2* + SX2)*(N*(SIN(A)*BIN*(M*(LKO*YK+SYLK)+LKO*PK*YK+SYPLK)+COS(A)*(M*( >LKO*XK+SXLK)+LKO*PK*XK+SXPLK))-(SPLK+M*SLK+LKO*(PK+M))*(SIN(A)*BIN ? *SY+COS(A)*SX))**2/(N*(SIN(A)**2*BIN**2*SY2+2*COS(A)*SIN(A)*BIN @ *SXY+COS(A)**2*SX2)-(SIN(A)*BIN*SY+COS(A)*SX)**2)**2+2*M*SPL2K+ 1 SP2L2K+M**2*SL2K+LKO**2*(PK+M)**2 RMSLK = SQRT(RMSLK/N) MINRMS(3+K) = RMSLK WRITE(TEXT,550) K,LK,LKOPT(K),RMSLK CALL STTPUT(TEXT,ISTAT) 50 CONTINUE ENDIF WRITE(TEXT,510) CALL STTPUT(TEXT,ISTAT) GOTO 2000 C ********************************************************************** 1000 CONTINUE C Computes RMS(A,M,BIN). There are too many parameters to create C a subroutine. Let's use the good old GOTO. RMSM = -N*((SPL+M*SL)*(SIN(A)**2*BIN**2*SY2+2*COS(A)*SIN(A)*BIN*S 1 XY+COS(A)**2*SX2)-(SIN(A)*BIN*SY+COS(A)*SX)*(SIN(A)*BIN*(SYPL+M 2 *SYL)+COS(A)*(SXPL+M*SXL)))**2/(N*(SIN(A)**2*BIN**2*SY2+2*COS(A 3 )*SIN(A)*BIN*SXY+COS(A)**2*SX2)-(SIN(A)*BIN*SY+COS(A)*SX)**2)** 4 2-2*(SIN(A)*BIN*SY+COS(A)*SX)*(N*(SIN(A)*BIN*(SYPL+M*SYL)+COS(A 5 )*(SXPL+M*SXL))-(SPL+M*SL)*(SIN(A)*BIN*SY+COS(A)*SX))*((SPL+M*S 6 L)*(SIN(A)**2*BIN**2*SY2+2*COS(A)*SIN(A)*BIN*SXY+COS(A)**2*SX2) 7 -(SIN(A)*BIN*SY+COS(A)*SX)*(SIN(A)*BIN*(SYPL+M*SYL)+COS(A)*(SXP 8 L+M*SXL)))/(N*(SIN(A)**2*BIN**2*SY2+2*COS(A)*SIN(A)*BIN*SXY+COS 9 (A)**2*SX2)-(SIN(A)*BIN*SY+COS(A)*SX)**2)**2-(SIN(A)**2*BIN**2* : SY2+2*COS(A)*SIN(A)*BIN*SXY+COS(A)**2*SX2)*(N*(SIN(A)*BIN*(SYPL ; +M*SYL)+COS(A)*(SXPL+M*SXL))-(SPL+M*SL)*(SIN(A)*BIN*SY+COS(A)*S < X))**2/(N*(SIN(A)**2*BIN**2*SY2+2*COS(A)*SIN(A)*BIN*SXY+COS(A)* = *2*SX2)-(SIN(A)*BIN*SY+COS(A)*SX)**2)**2+2*M*SPL2+SP2L2+M**2*SL > 2 RMSM = SQRT(RMSM/N) IF (MODE(1:1) .EQ. 'M') GOTO 1001 ! Absolute order number IF (MODE(1:1) .EQ. 'A') GOTO 1002 ! Angle IF (MODE(1:1) .EQ. 'B') GOTO 1003 ! Bin C ********************************************************************** 1900 CONTINUE ERROR = 1 2000 CONTINUE CALL STKWRI('OUTPUTI',ERROR,1,1,KUN,ISTAT) CALL STSEPI STOP END