C @(#)crfram.for 13.1.1.3 (ESO-DMD) 09/09/98 18:04:39 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 CRFRAM C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C program CRFRAM version 1.00 881122 C K. Banse ESO - Garching C 1.10 890512 1.20 900601 1.30 930110 C C.KEYWORDS C 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 new image file C DEFAULT/C/1/1 = Y, all image info given in command line C = N, use frame given in key P2 for that C P3/C/1/60 info frame, if DEFAULT = N C INPUTI/I/1/7 NAXIS + NPIX C INPUTR/R/1/12 START + STEP C P4/C/1/20 type of data function C = POLY, GAUSS, EXPO_DISK, RADIUS_LAW C = CIRCLE, ELLIPS C = NODATA if just a reference frame C = ASCII_FILE if input from a file C = TABLE_FILE if input from a table C P5/C/1/60 coefficients for data function C or name of data file if option ASCII above C or name of table file if option TABLE above C C.VERSIONS C 1.00 from version 3.60 (as of 880411) of CREFRAM C 1.10 throw away create/mask option C not necessary anymore (use copy/ii ...) C 1.20 add option to build 1dim frame from table C 1.30 add 3rd par. to CLNTAB call C 980706 C C-------------------------------------------------- C IMPLICIT NONE C INTEGER NAXIS,CUNLEN,MAPSIZ,IMSIZE INTEGER IAV,IDUM,METHOD,N,NCOEF,NO,STAT INTEGER INO,RNO,IBUF(7),NPIX(6) INTEGER IMNOX,EMPTY INTEGER UNIT(1),NLO,DFORM,MADRID(1) INTEGER*8 PNTR,PNTRX C REAL CUTS(4),COEFS(6),FMIN,FMAX C DOUBLE PRECISION START(6),STEP(6),DBUF(12),DDUM C CHARACTER CUNIT*112,IDENT*72,NAME*60,INFRAM*60 CHARACTER CBUF*60,DEFAUL*1 CHARACTER FUNCY*20,TABLE*60 C INCLUDE 'MID_INCLUDE:ST_DEF.INC' C COMMON /VMR/ MADRID C DATA CUNIT /'none given '/, CBUF /' '/ DATA NPIX /6*1/, START /6*0.D0/, STEP /6*1.D0/ DATA EMPTY /0/ DATA FMIN /+99999.99/ DATA FMAX /-99999.99/ DATA NAME /' '/, INFRAM /' '/ C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C get into MIDAS CALL STSPRO('CREFRAM') C C get name and type of new frame and DEFAULT flag CALL STKRDC('IN_A',1,1,60,IAV,NAME,UNIT,NLO,STAT) CALL STKRDC('DEFAULT',1,1,1,IAV,DEFAUL,UNIT,NLO,STAT) CALL UPCAS(DEFAUL,DEFAUL) C C image info from command line or reference frame IF (DEFAUL.NE.'Y') THEN CALL STKRDC('P3',1,1,60,IAV,INFRAM,UNIT,NLO,STAT) !get ref. frame C C read all relevant descriptors CALL STFOPN(INFRAM,D_OLD_FORMAT,0,F_IMA_TYPE,RNO,STAT) CALL STDRDI(RNO,'NAXIS',1,1,IAV,NAXIS,UNIT,NLO,STAT) IF (NAXIS.EQ.0) THEN EMPTY = 1 NAXIS = 1 ENDIF CALL STDRDI(RNO,'NPIX',1,NAXIS,IAV,NPIX,UNIT,NLO,STAT) CALL STDRDD(RNO,'START',1,NAXIS,IAV,START,UNIT,NLO,STAT) CALL STDRDD(RNO,'STEP',1,NAXIS,IAV,STEP,UNIT,NLO,STAT) CALL STDRDC(RNO,'IDENT',1,1,72,IAV,IDENT,UNIT,NLO,STAT) CALL STDRDC(RNO,'CUNIT',1,1,112,IAV,CUNIT,UNIT,NLO,STAT) ELSE C C get NAXIS, NPIX, START + STEP (max 6 dimensions) CALL STKRDI('INPUTI',1,7,IAV,IBUF,UNIT,NLO,STAT) CALL STKRDD('INPUTD',1,12,IAV,DBUF,UNIT,NLO,STAT) NAXIS = IBUF(1) C C set descriptor values IF ((NAXIS.LT.0).OR.(NAXIS.GT.6)) + CALL STETER(1,'invalid NAXIS ...') !invalid NAXIS ... C IF (NAXIS.EQ.0) THEN !want to build a real empty frame EMPTY = 1 NAXIS = 1 NPIX(1) = 1 START(1) = 0.0 STEP(1) = 1.0 ELSE DO 800, N=1,NAXIS NPIX(N) = IBUF(N+1) START(N) = DBUF(N) STEP(N) = DBUF(N+NAXIS) 800 CONTINUE ENDIF ENDIF C IMSIZE = 1 DO 900, N=1,NAXIS IMSIZE = IMSIZE * NPIX(N) 900 CONTINUE C C finally look for coefficients + function type CALL STKRDC('P4',1,1,20,IAV,FUNCY,UNIT,NLO,STAT) CALL UPCAS(FUNCY,FUNCY) CALL STKRDC('P5',1,1,60,IAV,CBUF,UNIT,NLO,STAT) C C be sure, to handle the NODATA option... DFORM = D_R4_FORMAT CUNLEN = (NAXIS+1)*16 C C create reference frame on disk + write descriptors IF (FUNCY(1:1).EQ.'N') THEN !NODATA option IMSIZE = 1 !frame with 1 pixel CALL STFCRE(NAME,DFORM,F_O_MODE,F_IMA_TYPE,IMSIZE,INO,STAT) CALL STDWRI(INO,'NAXIS',NAXIS,1,1,UNIT,STAT) CALL STDWRI(INO,'NPIX',NPIX,1,NAXIS,UNIT,STAT) CALL STDWRD(INO,'START',START,1,NAXIS,UNIT,STAT) CALL STDWRD(INO,'STEP',STEP,1,NAXIS,UNIT,STAT) IF (EMPTY.EQ.0) THEN IDENT = 'reference image (no data) ' ELSE IDENT = 'empty image (no data) ' ENDIF CALL STDWRC(INO,'IDENT',1,IDENT,1,72,UNIT,STAT) CALL STDWRC(INO,'CUNIT',1,CUNIT,1,CUNLEN,UNIT,STAT) GOTO 8800 !FMIN,FMAX are initialize to inhibit loading ELSE IDENT = 'artificial image ' !will be updated later CALL STKRDI('MONITPAR',20,1,IAV,MAPSIZ,UNIT,NLO,STAT) MAPSIZ = MAPSIZ * MAPSIZ IF ((NAXIS.GT.2) .AND. (MAPSIZ.LT.IMSIZE)) THEN CALL STFCRE(NAME,DFORM,F_O_MODE,F_IMA_TYPE,IMSIZE,INO,STAT) CALL STDWRI(INO,'NAXIS',NAXIS,1,1,UNIT,STAT) CALL STDWRI(INO,'NPIX',NPIX,1,NAXIS,UNIT,STAT) CALL STDWRD(INO,'START',START,1,NAXIS,UNIT,STAT) CALL STDWRD(INO,'STEP',STEP,1,NAXIS,UNIT,STAT) CALL STDWRC(INO,'IDENT',1,IDENT,1,72,UNIT,STAT) CALL STDWRC(INO,'CUNIT',1,CUNIT,1,CUNLEN,UNIT,STAT) PNTR = 0 !as reminder... ELSE CALL STIPUT(NAME,DFORM,F_O_MODE,F_IMA_TYPE, + NAXIS,NPIX,START,STEP,IDENT, + CUNIT(1:CUNLEN),PNTR,INO,STAT) ENDIF ENDIF C C branch according to function type IF (FUNCY(1:1).EQ.'P') THEN C C build polynomial of the form: p(x) = a + bx + cy + dxy + exx + fyy CALL GENCNV(CBUF,2,6,IDUM,COEFS,DDUM,NCOEF) IF (NCOEF.LE.0) THEN COEFS(1) = 0. NCOEF = 1 ENDIF C C and fill frame ... IF (PNTR.EQ.0) THEN !doti in pieces CALL STFCRE('fillwork',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE, + MAPSIZ,IMNOX,STAT) CALL STFMAP(IMNOX,F_X_MODE,1,MAPSIZ,IAV,PNTRX,STAT) CALL BIGFIL(MADRID(PNTRX),INO,IMSIZE,MAPSIZ,COEFS,FMIN,FMAX) CALL STFCLO(IMNOX,STAT) ELSE CALL IFILL(MADRID(PNTR),NAXIS,NPIX,START,STEP,COEFS, + NCOEF,FMIN,FMAX) ENDIF IDENT = 'image built from polygon ' GOTO 8800 C C input from an ASCII data file ELSE IF (FUNCY(1:1).EQ.'A') THEN IDENT = 'image created from ASCII file ' IDUM = 1 DO 3000 N=1,NAXIS IDUM = IDUM * NPIX(N) 3000 CONTINUE CALL DATFIL(CBUF,2,IDUM,MADRID(PNTR),MADRID(PNTR),1,FMIN,FMAX) IF (FMIN.GT.FMAX) CALL STTPUT + ('no valid data in file - pixels set to 0.0',STAT) GOTO 8800 C C input from a table with columns X:, Y: ELSE IF (FUNCY(1:1).EQ.'T') THEN IDENT = 'image created from Midas table ' CALL CLNTAB(CBUF,TABLE,0) CALL LINFOL(MADRID(PNTR),NPIX(1),START(1),STEP(1),TABLE) CALL MYMX(MADRID(PNTR),NPIX(1),CUTS) FMIN = CUTS(1) FMAX = CUTS(2) GOTO 8800 C ELSE IF (FUNCY(1:1).EQ.'G') THEN C C build Gaussian function IF ( (NAXIS.GE.3) .AND. (NPIX(3).GT.1) ) + CALL STETER(3,'not more than 2 dimensions supported yet') C NO = NAXIS * 2 !we expect NO numbers CALL GENCNV(CBUF,2,NO,IDUM,COEFS,DDUM,NCOEF) IF (NCOEF.LT.NO) THEN !default to center for mean, sigma = 1 COEFS(1) = START(1) + ((NPIX(1)-1)*STEP(1)*0.5) COEFS(2) = 3. * STEP(1) IF (NAXIS.EQ.2) THEN COEFS(3) = START(2) + ((NPIX(2)-1)*STEP(2)*0.5) COEFS(4) = 3. * STEP(2) ENDIF ENDIF C CALL GAUSS(MADRID(PNTR),NAXIS,NPIX,START,STEP,COEFS,NCOEF, + FMIN,FMAX) IF (NAXIS.EQ.1) THEN WRITE(IDENT,10000) COEFS(1),COEFS(2) ELSE WRITE(IDENT,10001) COEFS(1),COEFS(2),COEFS(3),COEFS(4) ENDIF GOTO 8800 ENDIF C CALL GENCNV(CBUF,2,4,IDUM,COEFS,DDUM,NCOEF) C C handle FFT filters IF (FUNCY(1:5).EQ.'FILT_') THEN IF (NCOEF.LT.3) COEFS(3) = 1. IF (NCOEF.LT.2) COEFS(2) = 1. IF (NCOEF.LT.1) COEFS(1) = ( STEP(1)+STEP(2) ) * 5 IF (COEFS(1).LT.10.E-12) + CALL STETER(2,'cut-off frequency must be > 0.0 ...') IF (FUNCY(6:7).EQ.'BL') METHOD = 5 IF (FUNCY(6:7).EQ.'EL') METHOD = 6 IF (FUNCY(6:7).EQ.'BH') METHOD = 7 IF (FUNCY(6:7).EQ.'EH') METHOD = 8 C C handle circle and ellips functions ELSE IF (FUNCY(1:1).EQ.'C') THEN METHOD = 3 IF (NCOEF.LT.3) COEFS(3) = 0. !outside IF (NCOEF.LT.2) COEFS(2) = 1. !inside IF (NCOEF.LT.1) COEFS(1) = !radius + (STEP(1)+STEP(2)) * (NPIX(1)+NPIX(2)) / 16. WRITE(IDENT,30001) COEFS(1),COEFS(2),COEFS(3) ELSE IF (FUNCY(1:2).EQ.'EL') THEN METHOD = 4 IF (NCOEF.LT.4) COEFS(4) = 0. !outside IF (NCOEF.LT.3) COEFS(3) = 1. !inside IF (NCOEF.LT.1) COEFS(1) = !a + (STEP(1)+STEP(2)) * (NPIX(1)+NPIX(2)) / 16. IF (NCOEF.LT.2) COEFS(2) = 0.6 * COEFS(1) !b WRITE(IDENT,40001) COEFS(1),COEFS(2),COEFS(3),COEFS(4) C C it remains exponential disk, radius law ELSE IF (NCOEF.LT.4) COEFS(4) = 0. IF (NCOEF.LT.3) COEFS(3) = 0. IF (NCOEF.LT.2) COEFS(2) = 4. IF (NCOEF.LT.1) COEFS(1) = 200. IF (FUNCY(1:1).EQ.'R') THEN METHOD = 2 WRITE(IDENT,50001) COEFS(1),COEFS(2),COEFS(3),COEFS(4) ELSE IF (FUNCY(1:2).NE.'EX') + CALL STTPUT + ('Invalid func_type, changed to `exponential disk''...', + STAT) METHOD = 1 !for exponential disk WRITE(IDENT,60001) COEFS(1),COEFS(2),COEFS(3),COEFS(4) ENDIF ENDIF CALL IMFUNC(METHOD,MADRID(PNTR),NAXIS,NPIX,START,STEP, + COEFS,FMIN,FMAX) C C update cuts and add descr. HISTORY 8800 CUTS(1) = 0. CUTS(2) = 0. CUTS(3) = FMIN CUTS(4) = FMAX CALL STDWRR(INO,'LHCUTS',CUTS,1,4,UNIT,STAT) CALL STDWRC(INO,'IDENT',1,IDENT,1,72,UNIT,STAT) CBUF(1:) = ' ' CALL DSCUPT(INO,INO,CBUF,STAT) C C check for empty frames IF (EMPTY.EQ.1) THEN NAXIS = 0 CALL STDWRI(INO,'NAXIS',NAXIS,1,1,UNIT,STAT) ENDIF C C that's it folks... CALL STSEPI C 10000 FORMAT('Gaussian, mean,sigma:',G12.5,',',G12.5) 10001 FORMAT('Gaussian, xm,xsig,ym,ysig:',3(G10.5,','),G10.5) 30001 FORMAT('Circle with coeffs =',2(G10.5,','),G10.5) 40001 FORMAT('Ellips with coeffs =',3(G10.5,','),G10.5) 50001 FORMAT('Radius Law, coeffs =',3(G10.5,','),G10.5) 60001 FORMAT('Exponential Disk, coeffs =',3(G10.5,','),G10.5) END SUBROUTINE BIGFIL(A,IMNO,IMSIZE,CHUNK,RC,RMIN,RMAX) C IMPLICIT NONE C INTEGER IMNO,IMSIZE,CHUNK INTEGER LOOPLM,RMAIND,STAT INTEGER N,K INTEGER MADRID(1) C REAL A(*),RC(*),RMIN,RMAX REAL RR C COMMON /VMR/ MADRID C LOOPLM = IMSIZE / CHUNK RMAIND = IMSIZE - (CHUNK*LOOPLM) K = 1 RR = RC(1) RMIN = RR RMAX = RR C DO 100, N=1,CHUNK A(N) = RR 100 CONTINUE C DO 1000, N=1,LOOPLM CALL STFPUT(IMNO,K,CHUNK,A,STAT) !data -> disk K = K + CHUNK 1000 CONTINUE C IF (RMAIND.NE.0) THEN CALL STFPUT(IMNO,K,RMAIND,A,STAT) ENDIF C RETURN END