C @(#)genran.for 13.1.1.1 (ESO-DMD) 06/02/98 18:17:30 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 GENRAN C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C program GENRAN version 1.00 850410 C F. MURTHAG STECF - Garching C M. ROSA STECF - Garching C K. Banse 1.20 860103 C K. Banse 1.50 881028 C C.KEYWORDS C random image C C.PURPOSE C create image filled with values taken from pseudo-random distributions C C.ALGORITHM C use NAGLIB C C.INPUT/OUTPUT C the following keys are used: C C OUT_A/C/1/60 result frame C IN_A/C/1/60 reference frame C ID/I/1/4 naxis,npix(*) C RD/D/1/6 start(*),step(*) C P4/C/1/20 function type C = GAUSSIAN, UNIFORM, EXPONENTIAL, C LOGNORMAL, BINOMIAL, POISSON, C CAUCHY C R/R/1/2 coefficients C I/I/1/1 seed C C.VERSION C C 1.50 move to FORTRAN 77 + new ST-Interfaces C C-------------------------------------------------- C IMPLICIT NONE C INTEGER IBUF(5),NPIX(3) INTEGER I,IAV,IROOT,NAXIS,NSIZE INTEGER*8 PNTR INTEGER IMNO,STAT INTEGER METHOD INTEGER UNI(1),NULO,MADRID(1) C REAL CUTS(4),COEFF(2) REAL FMIN,FMAX C DOUBLE PRECISION START(3),STEP(3),DBUF(6) C CHARACTER CUNIT*64,IDENTO*72,OUTFRM*60,REFIMA*60 CHARACTER CBUF*60 CHARACTER FNCTYP*20 C COMMON /VMR/ MADRID C INCLUDE 'MID_INCLUDE:ST_DEF.INC' C DATA IDENTO /'artificial image '/ DATA CUNIT /'none given '/ DATA CBUF /' '/ DATA NPIX /3*1/, START /3*0.D0/, STEP /3*1.D0/ DATA FMIN /+99999.99/ DATA FMAX /-99999.99/ C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C get into MIDAS C CALL STSPRO('GENRAN') C C get all keys C CALL STKRDC('OUT_A',1,1,60,IAV,OUTFRM,UNI,NULO,STAT) CALL STKRDI('I',1,1,IAV,IROOT,UNI,NULO,STAT) CALL STKRDR('R',1,2,IAV,COEFF,UNI,NULO,STAT) CALL STKRDC('P4',1,1,20,IAV,FNCTYP,UNI,NULO,STAT) CALL UPCAS(FNCTYP,FNCTYP) C C image info from command line or reference frame C CALL STKRDC('IN_A',1,1,60,IAV,REFIMA,UNI,NULO,STAT) IF (REFIMA(1:1).NE.'+') THEN CALL STFOPN(REFIMA,D_OLD_FORMAT,0,F_IMA_TYPE,IMNO,STAT) CALL STDRDI(IMNO,'NAXIS',1,1,IAV,NAXIS,UNI,NULO,STAT) CALL STDRDI(IMNO,'NPIX',1,NAXIS,IAV,NPIX,UNI,NULO,STAT) CALL STDRDD(IMNO,'START',1,NAXIS,IAV,START,UNI,NULO,STAT) CALL STDRDD(IMNO,'STEP',1,NAXIS,IAV,STEP,UNI,NULO,STAT) CALL STFCLO(IMNO,STAT) ELSE CALL STKRDI('ID',1,4,IAV,IBUF,UNI,NULO,STAT) CALL STKRDD('RD',1,6,IAV,DBUF,UNI,NULO,STAT) NAXIS = IBUF(1) IF ( (NAXIS.GE.1) .AND. (NAXIS.LE.3) ) THEN DO 400 I=1,NAXIS NPIX(I) = IBUF(I+1) START(I) = DBUF(I) STEP(I) = DBUF(NAXIS+I) 400 CONTINUE ELSE CALL STETER(1,'dimension not in [1,3]') ENDIF ENDIF C C create bulk data frame on disk + write descriptors C CALL STIPUT(OUTFRM,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXIS,NPIX,START,STEP,IDENTO, + CUNIT,PNTR,IMNO,STAT) C NSIZE=1 DO 600 I=1,NAXIS NSIZE = NSIZE*NPIX(I) 600 CONTINUE C METHOD = 1 IF (FNCTYP(1:1).EQ.'U') METHOD = 2 IF (FNCTYP(1:1).EQ.'E') METHOD = 3 IF (FNCTYP(1:1).EQ.'L') METHOD = 4 IF (FNCTYP(1:1).EQ.'B') METHOD = 5 IF (FNCTYP(1:1).EQ.'P') METHOD = 6 IF (FNCTYP(1:1).EQ.'C') METHOD = 7 C C generate random data C CALL RNDGEN(METHOD,COEFF,IROOT,MADRID(PNTR),NSIZE, + FMIN,FMAX) C C update cuts C CUTS (1) = FMIN CUTS (2) = FMAX CUTS (3) = FMIN CUTS (4) = FMAX CALL STDWRR(IMNO,'LHCUTS',CUTS,1,4,UNI,STAT) C C that's it folks... CALL STSEPI END SUBROUTINE RNDGEN(IFUNC,RCOF,ISEED,A,ISIZE,RMIN,RMAX) C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C Subroutine RANDGEN C Adapted by F. Murtagh from subroutine RANDOM in FUNCTIONS.FOR C by K. Banse, for *m*o*d*e*l*l*i*n*g*. 28.02.1985 C modified (cosmetics) M. Rosa 850410 C C.KEYWORDS C random number generation C C.PURPOSE C random number generation C C.ALGORITHM C C NAGLIB references: Subroutines and Functions C C S G05CBF(i) init repeatable random sequences C F G05DAF(a,b) real number between a,b from UNIFORM C F G05DBF(a) real number from NEG EXPONENTIAL with mean a C F G05DDF(a,b) real number from NORMAL with mean a and sigma b C F G05DEF(a,b) real number from LOG NORMAL, LN(x) mean a, sigma b C F G05DFF(a,,b) real number from CAUCHY, median a, semi-interq. range b C S G05ECF(t,,,) ref. vector for POISSON with mean t(real) C S G05EDF(n,p,,,) ref. vector for BINOMIAL with n trials with prob. p C F G05EYF(,) return random integer from discrete distributions G05ExF C C.INPUT/OUTPUT C C-------------------------------------------------- C IMPLICIT NONE C INTEGER IFUNC,ISEED,ISIZE INTEGER I,IFAIL,NTRY,NR,NX,STAT INTEGER G05EYF C REAL A(*),RCOF(*),RMIN,RMAX REAL PIX C DOUBLE PRECISION RC1,RC2 DOUBLE PRECISION T,REF(2000) DOUBLE PRECISION G05DAF,G05DBF,G05DDF,G05DEF,G05DFF C CHARACTER OUTPUT*80 C RC1 = RCOF(1) RC2 = RCOF(2) C C init the random number generator CALL G05CBF(ISEED) C C branch according to method GOTO (100,200,300,400,500,600,700),IFUNC C C Gaussian distribution C 100 DO 130 I=1,ISIZE PIX = G05DDF(RC1,RC2) IF (RMIN.GT.PIX) RMIN = PIX IF (RMAX.LT.PIX) RMAX = PIX A(I) = PIX 130 CONTINUE RETURN C C Uniform distribution C 200 DO 230 I=1,ISIZE PIX = G05DAF(RC1,RC2) IF (RMIN.GT.PIX) RMIN = PIX IF (RMAX.LT.PIX) RMAX = PIX A(I) = PIX 230 CONTINUE RETURN C C (Negative) exponential distribution C 300 DO 330 I=1,ISIZE PIX = G05DBF(RC1) IF (RMIN.GT.PIX) RMIN = PIX IF (RMAX.LT.PIX) RMAX = PIX A(I) = PIX 330 CONTINUE RETURN C C Lognormal distribution C 400 DO 430 I=1,ISIZE PIX = G05DEF(RC1,RC2) IF (RMIN.GT.PIX) RMIN = PIX IF (RMAX.LT.PIX) RMAX = PIX A(I) = PIX 430 CONTINUE RETURN C C Binomial distribution C 500 IFAIL = 0 NTRY = 100 NR = 2000 CALL G05EDF(NTRY,RC1,REF,NR,IFAIL) !compute reference vector IF (IFAIL.NE.0) THEN WRITE(OUTPUT,10000) IFAIL CALL STTPUT(OUTPUT,STAT) RETURN ENDIF DO 530 I=1,ISIZE NX = G05EYF(REF,NR) PIX = FLOAT(NX) IF (RMIN.GT.PIX) RMIN = PIX IF (RMAX.LT.PIX) RMAX = PIX A(I) = PIX 530 CONTINUE RETURN C C Poisson distribution C 600 IFAIL = 0 T = RC1 IF (T.LT.0.) T = 0. IF (T.GT.10000.) T = 10000. C T = MIN(MAX(0.,RC1),10000.) !force T into [0,10000] NR = 2000 CALL G05ECF(T,REF,NR,IFAIL) !compute reference vector IF (IFAIL.NE.0) THEN WRITE(OUTPUT,10000) IFAIL CALL STTPUT(OUTPUT,STAT) RETURN ENDIF DO 630 I=1,ISIZE NX = G05EYF(REF,NR) PIX = FLOAT(NX) IF (RMIN.GT.PIX) RMIN = PIX IF (RMAX.LT.PIX) RMAX = PIX A(I) = PIX 630 CONTINUE RETURN C C Cauchy distribution C 700 DO 730 I=1,ISIZE PIX = G05DFF(RC1,RC2) IF (RMIN.GT.PIX) RMIN = PIX IF (RMAX.LT.PIX) RMAX = PIX A(I) = PIX 730 CONTINUE RETURN C 10000 FORMAT(' error return from NAG lib. = ',I1) END