C @(#)sdp.for	19.1 (ESO-DMD) 02/25/03 13:24:48
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
      PROGRAM SDP
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  program SDP			version 1.00	910110
C  K. Banse				ESO - Garching
C    1.20	920430
C
C.KEYWORDS
C  dynamic systems, bulk data frame
C
C.PURPOSE
C  create a MIDAS image
C
C.ALGORITHM
C  use MIDAS interfaces
C
C.INPUT/OUTPUT
C  the following keywords are used:
C
C IN_A/C/1/60                  	name of input 1-dim image
C OUT_A/C/1/60                 	name of output 2-dim image
C
C.VERSIONS
C
C--------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER      NAXIS
      INTEGER*8    PNTRA,PNTRB,PNTRW
      INTEGER      IAV,NO,STAT,IBUF(10)
      INTEGER      INOA,INOB,INOW,NPIXA,NPIX(2)
      INTEGER      NTOP,NLO,UNIT(1),MADRID(1)
C
      DOUBLE PRECISION START(2),STEP(2)
C
      CHARACTER    CUNIT*48,IDENT*72,INFRA*60,OUTFRA*60
      CHARACTER    CBUF*60
C 
      REAL         CUTS(4)
C
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
C
      COMMON   /VMR/ MADRID
C
      DATA      CUNIT /'none given '/, CBUF /' '/
      DATA      NPIX      /2*1/,       START /2*0.D0/, STEP /2*1.D0/
C
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
C  get into MIDAS
      CALL STSPRO('SDP')
C
C  get name of input, result frame
      CALL STKRDC('IN_A',1,1,60,IAV,INFRA,UNIT,NLO,STAT)
      CALL STKRDC('OUT_A',1,1,60,IAV,OUTFRA,UNIT,NLO,STAT)
      CALL STKRDI('INPUTI',1,10,IAV,IBUF,UNIT,NLO,STAT)
      NTOP = IBUF(1)
C
C  get input frame + create working buffer and output graph
      CALL STIGET(INFRA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     +            1,NAXIS,NPIX,START,STEP,IDENT,
     +            CUNIT,PNTRA,INOA,STAT)
      NAXIS = 1                                   !force to 1-dim frame...
      CALL STIPUT('dummy',D_I4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +            NAXIS,NPIX,START,STEP,IDENT,
     +            CUNIT,PNTRW,INOW,STAT)
      NPIXA = NPIX(1)
      CALL STDRDR(INOA,'LHCUTS',1,4,IAV,CUTS,UNIT,NLO,STAT)
      IF (CUTS(1).GE.CUTS(2)) THEN
         NO = 3
      ELSE
         NO = 1
      ENDIF
C 
      NPIX(1) = 2*NTOP + 1
      NPIX(2) = NPIX(1)
      NAXIS = 2
      START(1) = -NTOP
      START(2) = START(1)
      STEP(1) = 1.
      STEP(2) = STEP(1)
      CALL STIPUT(OUTFRA,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +            NAXIS,NPIX,START,STEP,IDENT,
     +            CUNIT,PNTRB,INOB,STAT)
C
C map dat into interval [0,NTOP]
      CALL MYSCAL(MADRID(PNTRA),MADRID(PNTRW),NPIXA,NTOP,CUTS(NO))
C
C and do the symmetrical dot pattern
      CALL DOSDP(MADRID(PNTRW),NPIXA,MADRID(PNTRB),NPIX,IBUF)
C 
C  that's it folks...
      CALL STSEPI
      END

      SUBROUTINE MYSCAL(A,IB,NDIM,NTOP,FCUTS)
C 
      IMPLICIT NONE
C 
      INTEGER   IB(1),NDIM,NTOP
      INTEGER   N
C 
      REAL      A(1),FCUTS(2)
      REAL      AMIN,AMAX,FACT
C 
      AMIN = FCUTS(1)
      AMAX = FCUTS(2)
C 
C  now scale
      FACT = NTOP / (AMAX-AMIN)
      DO 2000 N=1,NDIM
         IF (A(N).LE.AMIN) THEN
            IB(N) = 0
         ELSE IF (A(N).GE.AMAX) THEN
            IB(N) = NTOP
         ELSE
            IB(N) = NINT((A(N)-AMIN)*FACT)
         ENDIF
2000  CONTINUE
C 
      RETURN
      END

      SUBROUTINE DOSDP(IB,NDIM,C,NPIX,KBUF)
C 
      IMPLICIT NONE
C 
      INTEGER  IB(*),NDIM,NPIX(2),KBUF(*)
      INTEGER  N,NN,KR,JB,NHALF,OFF
      INTEGER  LAG,IANG,IX,IY
C 
      REAL  C(1)
      REAL  PI,RF,ANGLE,COLOR,CO,SI
C 
      DATA     PI /3.141593/                         !Pi
C 
      DO 200 N=1,NPIX(1)*NPIX(2)
         C(N) = 0.0
200   CONTINUE
C 
      IANG = 60
      LAG = KBUF(2)
      RF = PI / 180.
      NHALF = KBUF(1) + 1
C 
      DO 1000 N=1,NDIM-LAG
         JB = IB(N+LAG)
         KR = IB(N)
         IF (JB.GE.KR) THEN               !IB(n+lag) >= IB(n)
            COLOR = KBUF(3)
         ELSE
            COLOR = KBUF(4)
         ENDIF
C
         DO 500 NN=1,360,IANG
            IF (KR.EQ.0) THEN
               OFF = (NHALF-1)*NPIX(1) + NHALF
               IF (KBUF(5).EQ.1) THEN
                  C(OFF) = C(OFF) + COLOR 
               ELSE
                  C(OFF) = COLOR
               ENDIF
            ELSE
               ANGLE = NN + JB
               ANGLE = ANGLE * RF
               CO = COS(ANGLE)
               SI = SIN(ANGLE)
               IX = NINT(KR * CO)
               IY = NINT(KR * SI)
               IX = IX + NHALF                  !center of array = 0
               IY = IY + NHALF 
               OFF = (IY-1)*NPIX(1) + IX
               IF (KBUF(5).EQ.1) THEN
                  C(OFF) = C(OFF) + COLOR
               ELSE
                  C(OFF) = COLOR
               ENDIF
C 
               ANGLE = NN - JB
               ANGLE = ANGLE * RF
               CO = COS(ANGLE)
               SI = SIN(ANGLE)
               IX = NINT(KR * CO)
               IY = NINT(KR * SI)
               IX = IX + NHALF                  !center of array = 0
               IY = IY + NHALF 
               OFF = (IY-1)*NPIX(1) + IX
               IF (KBUF(5).EQ.1) THEN
                  C(OFF) = C(OFF) + COLOR
               ELSE
                  C(OFF) = COLOR
               ENDIF
            ENDIF
500      CONTINUE
1000  CONTINUE
C 
      RETURN
      END
