C @(#)dcalc.for 17.1.1.1 (ESO-DMD) 01/25/02 17:40:33 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 Massachusetts Ave, Cambridge, C MA 02139, USA. C C Correspondence 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 SUBROUTINE DOPCC(OPCODE,CONSTA,CONSTB,CONSTC) C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subroutine OPCC version 2.00 880202 C K. Banse ESO - Garching C source code on file CALC.F C 2.10 890918 C C.KEYWORDS C bulk data frames, arithmetic operations C C.PURPOSE C perform arithmetic on arrays C C.ALGORITHM C do it according to opcode ( ", +, -, *, / or =) C " means ** (to the power of) C special condition handler is used for division by "zero" C C.INPUT/OUTPUT C call as OPCC(OPCODE,CONSTA,CONSTB,CONSTC) C C input par: C OPCODE: char.exp. operation code in postfix notation, e.g. "FF+" C A: R*8 array 1. input array C B: R*8 array 2. input array C SCALR: R*8 scalar value C CONSTA: R*8 1. input scalar C CONSTB: R*8 2. input scalar C NDIM: I*4 size of arrays involved C C output par: C C: R*8 array result array C CONSTC: R*8 result scalar C C-------------------------------------------------- C IMPLICIT NONE C INTEGER NCOUNT,TSTFLG C DOUBLE PRECISION CONSTA,CONSTB,CONSTC DOUBLE PRECISION USRNUL,DEPSP,DEPSN C CHARACTER*(*) OPCODE C COMMON /NULCOM/ NCOUNT,TSTFLG,USRNUL,DEPSP,DEPSN C C branch according to opcode C C subtract 2. scalar from 1. scalar IF (OPCODE(3:3).EQ.'-') THEN CONSTC = CONSTA - CONSTB C C add scalars ELSE IF (OPCODE(3:3).EQ.'+') THEN CONSTC = CONSTA + CONSTB C C divide 1. scalar by 2. scalar ELSE IF (OPCODE(3:3).EQ.'/') THEN IF ((CONSTB.LT.DEPSP) .AND. + (CONSTB.GT.DEPSN)) THEN !always check divisor... CONSTC = USRNUL NCOUNT = NCOUNT + 1 ELSE CONSTC = CONSTA / CONSTB ENDIF C C multiply scalars ELSE IF (OPCODE(3:3).EQ.'*') THEN CONSTC = CONSTA * CONSTB C C compute 1. scalar ** 2. scalar ELSE IF (OPCODE(3:3).EQ.'`') THEN CONSTC = CONSTA ** CONSTB ENDIF C RETURN END SUBROUTINE DOPCF(A,CVAL,NDIM) C IMPLICIT NONE C INTEGER N,NDIM C DOUBLE PRECISION A(*),CVAL C DO 1000, N=1,NDIM A(N) = CVAL 1000 CONTINUE C RETURN END SUBROUTINE DF1C(CFUNC,AVAL,CVAL) C IMPLICIT NONE C INTEGER NCOUNT,TSTFLG C DOUBLE PRECISION AVAL,CVAL DOUBLE PRECISION TEMP,USRNUL,DEPSP,DEPSN,PI C CHARACTER*5 CFUNC C COMMON /NULCOM/ NCOUNT,TSTFLG,USRNUL,DEPSP,DEPSN C DATA PI /3.1415926535/ !Pi C C branch according to function C C SQRT IF (CFUNC.EQ.'SQRT ') THEN IF (AVAL.LT.0.D0) THEN CVAL = DBLE(USRNUL) NCOUNT = NCOUNT + 1 ELSE CVAL = DSQRT(AVAL) ENDIF C C LN ELSE IF (CFUNC.EQ.'LN ') THEN IF (AVAL.LE.0.D0) THEN CVAL = USRNUL NCOUNT = NCOUNT + 1 ELSE CVAL = DLOG(AVAL) ENDIF C C LOG10 or LOG ELSE IF (CFUNC(1:3).EQ.'LOG') THEN IF (AVAL.LE.0.D0) THEN CVAL = USRNUL NCOUNT = NCOUNT + 1 ELSE CVAL = DLOG10(AVAL) ENDIF C C EXP ELSE IF (CFUNC.EQ.'EXP ') THEN CVAL = DEXP(AVAL) C C EXP10 ELSE IF (CFUNC.EQ.'EXP10') THEN CVAL = 10.D0 ** AVAL C C SIN ELSE IF (CFUNC.EQ.'SIN ') THEN TEMP = (AVAL / 180.D0) * PI !degrees -> radians CVAL = DSIN(TEMP) C C COS ELSE IF (CFUNC.EQ.'COS ') THEN TEMP = (AVAL / 180.D0) * PI !degrees -> radians CVAL = DCOS(TEMP) C C TAN ELSE IF (CFUNC.EQ.'TAN ') THEN TEMP = (AVAL / 180.D0) * PI !degrees -> radians CVAL = DTAN(TEMP) C C ASIN ELSE IF (CFUNC.EQ.'ASIN ') THEN TEMP = DASIN(AVAL) CVAL = (TEMP * 180.D0) / PI !radians -> degrees C C ACOS ELSE IF (CFUNC.EQ.'ACOS ') THEN TEMP = DACOS(AVAL) CVAL = (TEMP * 180.D0) / PI !radians -> degrees C C ATAN ELSE IF (CFUNC.EQ.'ATAN ') THEN TEMP = DATAN(AVAL) CVAL = (TEMP * 180.D0) / PI !radians -> degrees C C NINT ELSE IF (CFUNC.EQ.'INT ') THEN CVAL = DBLE(NINT(AVAL)) C C ABS ELSE IF (CFUNC.EQ.'ABS ') THEN CVAL = DABS(AVAL) ENDIF C RETURN END SUBROUTINE DF2CC(CFUNC,AVAL,CVAL) C IMPLICIT NONE C DOUBLE PRECISION AVAL(2),CVAL,PI DOUBLE PRECISION TEMP C CHARACTER*5 CFUNC C DATA PI /3.1415926535/ !Pi C C branch according to function C C function ATAN2 IF (CFUNC.EQ.'ATAN2') THEN TEMP = DATAN2(AVAL(1),AVAL(2)) CVAL = (TEMP * 180.D0) / PI !radians -> degrees C C get minimum ELSE IF (CFUNC.EQ.'MIN ') THEN CVAL = MIN(AVAL(1),AVAL(2)) C C get maximum ELSE IF (CFUNC.EQ.'MAX ') THEN CVAL = MAX(AVAL(1),AVAL(2)) C C get remainder ELSE IF (CFUNC.EQ.'MOD ') THEN CVAL = DMOD(AVAL(1),AVAL(2)) ENDIF C RETURN END SUBROUTINE DOPFW(OPCODE,A,B,C,APIX,BPIX,CPIX, + NPIXA,NPIXB,NPIXC) C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subroutine OPFFW version 2.00 880202 C FN2FFW 2.00 880202 C K. Banse ESO - Garching C 2.10 890915 C C.KEYWORDS C bulk data frames, arithmetic operations C C.PURPOSE C perform arithmetic on arrays C C.ALGORITHM C do it according to opcode ( `, +, -, *, / or =) C ` means ** (to the power of) C special condition handler is used for division by "zero" C C.INPUT/OUTPUT C call as OPFFW(OPCODE,A,B,C,APIX,BPIX,CPIX,NPIXA,NPIXB,NPIXC) C C input par: C OPCODE: char.exp. operation code in postfix notation, e.g. "FF+" C A: R*8 array 1. input array C B: R*8 array 2. input array C APIX: I*4 array start + end pixels for A C BPIX: I*4 array start + end pixels for B C CPIX: I*4 array start + end pixels for C C NPIXA: I*4 no. of pixels per line in A C NPIXB: I*4 no. of pixels per line in B C NPIXC: I*4 no. of pixels per line in C C C output par: C C: R*8 array result array C C-------------------------------------------------- C IMPLICIT NONE C INTEGER N,NN,NNN,NCOUNT,TSTFLG INTEGER NX,NY,NZ,OFFA,OFFB,OFFC INTEGER XOFFA,XOFFB,XOFFC INTEGER YOFFA,YOFFB,YOFFC INTEGER ZOFFA,ZOFFB,ZOFFC INTEGER APIX(3,2),BPIX(3,2),CPIX(3,2) INTEGER NPIXA(*),NPIXB(*),NPIXC(*) C CHARACTER*(*) OPCODE C DOUBLE PRECISION A(*),B(*),C(*) DOUBLE PRECISION USRNUL,DEPSP,DEPSN C COMMON /NULCOM/ NCOUNT,TSTFLG,USRNUL,DEPSP,DEPSN C C init NX = APIX(1,2) - APIX(1,1) + 1 NY = APIX(2,2) - APIX(2,1) + 1 NZ = APIX(3,2) - APIX(3,1) + 1 C XOFFA = APIX(1,1) - 1 XOFFC = CPIX(1,1) - 1 YOFFA = (APIX(2,1)-1)*NPIXA(1) YOFFC = (CPIX(2,1)-1)*NPIXC(1) ZOFFA = (APIX(3,1)-1)*NPIXA(1)*NPIXA(2) ZOFFC = (CPIX(3,1)-1)*NPIXC(1)*NPIXC(2) C C copy 1. frame into result frame IF (OPCODE(3:3).EQ.'=') THEN DO 5500, NNN=1,NZ OFFA = ZOFFA + YOFFA + XOFFA OFFC = ZOFFC + YOFFC + XOFFC C DO 5300, NN = 1,NY DO 5200, N=1,NX C(OFFC+N) = A(OFFA+N) 5200 CONTINUE OFFA = OFFA + NPIXA(1) OFFC = OFFC + NPIXC(1) 5300 CONTINUE C ZOFFA = ZOFFA + NPIXA(1)*NPIXA(2) ZOFFC = ZOFFC + NPIXC(1)*NPIXC(2) 5500 CONTINUE RETURN ENDIF C C XOFFB = BPIX(1,1) - 1 !all other opcodes have 3 operands YOFFB = (BPIX(2,1)-1)*NPIXB(1) ZOFFB = (BPIX(3,1)-1)*NPIXB(1)*NPIXB(2) C C branch according to opcode C C add frames IF (OPCODE(3:3).EQ.'+') THEN DO 1500, NNN=1,NZ OFFA = ZOFFA + YOFFA + XOFFA OFFB = ZOFFB + YOFFB + XOFFB OFFC = ZOFFC + YOFFC + XOFFC C DO 1300, NN = 1,NY DO 1100, N=1,NX C(OFFC+N) = A(OFFA+N) + B(OFFB+N) 1100 CONTINUE OFFA = OFFA + NPIXA(1) OFFB = OFFB + NPIXB(1) OFFC = OFFC + NPIXC(1) 1300 CONTINUE C ZOFFA = ZOFFA + NPIXA(1)*NPIXA(2) ZOFFB = ZOFFB + NPIXB(1)*NPIXB(2) ZOFFC = ZOFFC + NPIXC(1)*NPIXC(2) 1500 CONTINUE C C subtract 2. frame from 1. frame ELSE IF (OPCODE(3:3).EQ.'-') THEN DO 2500, NNN=1,NZ OFFA = ZOFFA + YOFFA + XOFFA OFFB = ZOFFB + YOFFB + XOFFB OFFC = ZOFFC + YOFFC + XOFFC C DO 2300, NN = 1,NY DO 2200, N=1,NX C(OFFC+N) = A(OFFA+N) - B(OFFB+N) 2200 CONTINUE OFFA = OFFA + NPIXA(1) OFFB = OFFB + NPIXB(1) OFFC = OFFC + NPIXC(1) 2300 CONTINUE C ZOFFA = ZOFFA + NPIXA(1)*NPIXA(2) ZOFFB = ZOFFB + NPIXB(1)*NPIXB(2) ZOFFC = ZOFFC + NPIXC(1)*NPIXC(2) 2500 CONTINUE C C multiply frames ELSE IF (OPCODE(3:3).EQ.'*') THEN DO 3500, NNN=1,NZ OFFA = ZOFFA + YOFFA + XOFFA OFFB = ZOFFB + YOFFB + XOFFB OFFC = ZOFFC + YOFFC + XOFFC C DO 3300, NN = 1,NY DO 3200, N=1,NX C(OFFC+N) = A(OFFA+N) * B(OFFB+N) 3200 CONTINUE OFFA = OFFA + NPIXA(1) OFFB = OFFB + NPIXB(1) OFFC = OFFC + NPIXC(1) 3300 CONTINUE C ZOFFA = ZOFFA + NPIXA(1)*NPIXA(2) ZOFFB = ZOFFB + NPIXB(1)*NPIXB(2) ZOFFC = ZOFFC + NPIXC(1)*NPIXC(2) 3500 CONTINUE C C divide 1. frame by 2. frame ELSE IF (OPCODE(3:3).EQ.'/') THEN C IF (TSTFLG.EQ.1) THEN DO 4500, NNN=1,NZ OFFA = ZOFFA + YOFFA + XOFFA OFFB = ZOFFB + YOFFB + XOFFB OFFC = ZOFFC + YOFFC + XOFFC C DO 4300, NN = 1,NY DO 4200, N=1,NX IF ((B(OFFB+N) .LT. DEPSP) .AND. + (B(OFFB+N) .GT. DEPSN)) THEN C(OFFC+N) = USRNUL NCOUNT = NCOUNT + 1 ELSE C(OFFC+N) = A(OFFA+N) / B(OFFB+N) ENDIF 4200 CONTINUE OFFA = OFFA + NPIXA(1) OFFB = OFFB + NPIXB(1) OFFC = OFFC + NPIXC(1) 4300 CONTINUE C ZOFFA = ZOFFA + NPIXA(1)*NPIXA(2) ZOFFB = ZOFFB + NPIXB(1)*NPIXB(2) ZOFFC = ZOFFC + NPIXC(1)*NPIXC(2) 4500 CONTINUE ELSE C DO 4800, NNN=1,NZ OFFA = ZOFFA + YOFFA + XOFFA OFFB = ZOFFB + YOFFB + XOFFB OFFC = ZOFFC + YOFFC + XOFFC C DO 4700, NN = 1,NY DO 4600, N=1,NX C(OFFC+N) = A(OFFA+N) / B(OFFB+N) 4600 CONTINUE OFFA = OFFA + NPIXA(1) OFFB = OFFB + NPIXB(1) OFFC = OFFC + NPIXC(1) 4700 CONTINUE C ZOFFA = ZOFFA + NPIXA(1)*NPIXA(2) ZOFFB = ZOFFB + NPIXB(1)*NPIXB(2) ZOFFC = ZOFFC + NPIXC(1)*NPIXC(2) 4800 CONTINUE ENDIF C C compute 1. frame ** 2. frame ELSE IF (OPCODE(3:3).EQ.'`') THEN DO 6500, NNN=1,NZ OFFA = ZOFFA + YOFFA + XOFFA OFFB = ZOFFB + YOFFB + XOFFB OFFC = ZOFFC + YOFFC + XOFFC C DO 6300, NN = 1,NY DO 6200, N=1,NX C(OFFC+N) = A(OFFA+N) ** B(OFFB+N) 6200 CONTINUE OFFA = OFFA + NPIXA(1) OFFB = OFFB + NPIXB(1) OFFC = OFFC + NPIXC(1) 6300 CONTINUE C ZOFFA = ZOFFA + NPIXA(1)*NPIXA(2) ZOFFB = ZOFFB + NPIXB(1)*NPIXB(2) ZOFFC = ZOFFC + NPIXC(1)*NPIXC(2) 6500 CONTINUE ENDIF C RETURN END SUBROUTINE DF2FFW(CFUNC,A,B,C,APIX,BPIX,CPIX, + NPIXA,NPIXB,NPIXC) C IMPLICIT NONE C INTEGER N,NN,NNN,NCOUNT,TSTFLG INTEGER NX,NY,NZ,OFFA,OFFB,OFFC INTEGER XOFFA,XOFFB,XOFFC INTEGER YOFFA,YOFFB,YOFFC INTEGER ZOFFA,ZOFFB,ZOFFC INTEGER APIX(3,2),BPIX(3,2),CPIX(3,2) INTEGER NPIXA(*),NPIXB(*),NPIXC(*) C CHARACTER*5 CFUNC C DOUBLE PRECISION A(*),B(*),C(*),FACT,RA,RB DOUBLE PRECISION USRNUL,DEPSP,DEPSN C COMMON /NULCOM/ NCOUNT,TSTFLG,USRNUL,DEPSP,DEPSN C DATA FACT /0.0174532925/ C C init NX = APIX(1,2) - APIX(1,1) + 1 NY = APIX(2,2) - APIX(2,1) + 1 NZ = APIX(3,2) - APIX(3,1) + 1 C XOFFA = APIX(1,1) - 1 XOFFB = BPIX(1,1) - 1 XOFFC = CPIX(1,1) - 1 YOFFA = (APIX(2,1)-1)*NPIXA(1) YOFFB = (BPIX(2,1)-1)*NPIXB(1) YOFFC = (CPIX(2,1)-1)*NPIXC(1) ZOFFA = (APIX(3,1)-1)*NPIXA(1)*NPIXA(2) ZOFFB = (BPIX(3,1)-1)*NPIXB(1)*NPIXB(2) ZOFFC = (CPIX(3,1)-1)*NPIXC(1)*NPIXC(2) C C branch according to function C C function ATAN2 IF (CFUNC.EQ.'ATAN2') THEN DO 1500, NNN=1,NZ OFFA = ZOFFA + YOFFA + XOFFA OFFB = ZOFFB + YOFFB + XOFFB OFFC = ZOFFC + YOFFC + XOFFC C DO 1300, NN = 1,NY DO 1200, N=1,NX RA = A(OFFA+N) * FACT RB = B(OFFB+N) * FACT C(OFFC+N) = DATAN2(RA,RB) 1200 CONTINUE OFFA = OFFA + NPIXA(1) OFFB = OFFB + NPIXB(1) OFFC = OFFC + NPIXC(1) 1300 CONTINUE C ZOFFA = ZOFFA + NPIXA(1)*NPIXA(2) ZOFFB = ZOFFB + NPIXB(1)*NPIXB(2) ZOFFC = ZOFFC + NPIXC(1)*NPIXC(2) 1500 CONTINUE C C minimum of 2 frames ELSE IF (CFUNC.EQ.'MIN ') THEN DO 2500, NNN=1,NZ OFFA = ZOFFA + YOFFA + XOFFA OFFB = ZOFFB + YOFFB + XOFFB OFFC = ZOFFC + YOFFC + XOFFC C DO 2300, NN = 1,NY DO 2200, N=1,NX C(OFFC+N) = MIN(A(OFFA+N),B(OFFB+N)) 2200 CONTINUE OFFA = OFFA + NPIXA(1) OFFB = OFFB + NPIXB(1) OFFC = OFFC + NPIXC(1) 2300 CONTINUE C ZOFFA = ZOFFA + NPIXA(1)*NPIXA(2) ZOFFB = ZOFFB + NPIXB(1)*NPIXB(2) ZOFFC = ZOFFC + NPIXC(1)*NPIXC(2) 2500 CONTINUE C C maximum of 2 frames ELSE IF (CFUNC.EQ.'MAX ') THEN DO 3500, NNN=1,NZ OFFA = ZOFFA + YOFFA + XOFFA OFFB = ZOFFB + YOFFB + XOFFB OFFC = ZOFFC + YOFFC + XOFFC C DO 3300, NN = 1,NY DO 3200, N=1,NX C(OFFC+N) = MAX(A(OFFA+N),B(OFFB+N)) 3200 CONTINUE OFFA = OFFA + NPIXA(1) OFFB = OFFB + NPIXB(1) OFFC = OFFC + NPIXC(1) 3300 CONTINUE C ZOFFA = ZOFFA + NPIXA(1)*NPIXA(2) ZOFFB = ZOFFB + NPIXB(1)*NPIXB(2) ZOFFC = ZOFFC + NPIXC(1)*NPIXC(2) 3500 CONTINUE C C remainder of 1.frame/2. frames ELSE IF (CFUNC.EQ.'MOD ') THEN DO 4500, NNN=1,NZ OFFA = ZOFFA + YOFFA + XOFFA OFFB = ZOFFB + YOFFB + XOFFB OFFC = ZOFFC + YOFFC + XOFFC C DO 4300, NN = 1,NY DO 4200, N=1,NX C(OFFC+N) = DMOD(A(OFFA+N),B(OFFB+N)) 4200 CONTINUE OFFA = OFFA + NPIXA(1) OFFB = OFFB + NPIXB(1) OFFC = OFFC + NPIXC(1) 4300 CONTINUE C ZOFFA = ZOFFA + NPIXA(1)*NPIXA(2) ZOFFB = ZOFFB + NPIXB(1)*NPIXB(2) ZOFFC = ZOFFC + NPIXC(1)*NPIXC(2) 4500 CONTINUE ENDIF C RETURN END SUBROUTINE DDSPNL(NCC) C INTEGER NCC INTEGER NCOUNT,TSTFLG,STAT C CHARACTER OUTPUT*80 C DOUBLE PRECISION USRNUL,DEPSP,DEPSN C COMMON /NULCOM/ NCOUNT,TSTFLG,USRNUL,DEPSP,DEPSN C IF (NCC.GT.1) THEN WRITE(OUTPUT,10000) NCC,USRNUL ELSE WRITE(OUTPUT,10001) USRNUL ENDIF CALL STTPUT(OUTPUT,STAT) RETURN C C format C 10000 FORMAT(I7,' undefined pixels ... set to "null value" = ',G15.7) 10001 FORMAT('1 undefined pixel ... set to "null value" = ',G15.7) C END