cc lsfd.f program lsfd implicit real*8 (a-h,o-z) cc ================================================ cc Set parameters: cc ~~~~~~~~~~~~~~~ cc nsim : number of simulated light curves cc (in Monte Carlo simulation) parameter(nsim = 100) cc asig : threshold value for sigma-clipping parameter(asig = 3.0d0) cc ================================================ cc Explicit declarations: cc ~~~~~~~~~~~~~~~~~~~~~~ character * ( 80 ) arg(100) character filets * (80) character fileprw * (80) integer iargc integer numarg real*8 w(100) real*8 t(10000) real*8 x(10000) real*8 xw(10000) real*8 amp(1000,100), phs(1000,100),sfreq(1000) real*8 aa(100), pp(100) real*8 a(100),b(100) real*8 ampavr(100) real*8 ampdsp(100) real*8 phsavr(100) real*8 phsdsp(100) cc ================================================ cc Set basic logical parameters: cc ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ cc iwhite = 1 --> Write out the prewhitened cc time series into 'lsfd.prw' iwhite = 0 cc ierror = 1 --> Calculate the amplitude and cc phase errors by Monte Carlo simulation ierror = 1 cc ifitf = 1 --> Fit the frequency. cc ifitf = 1 p7 = 8.0d0*datan(1.0d0) cc cc dfra = accuracy of the frequency determination in terms cc of the line width: dfra=df/(1/T) [USED IF IFITF=1] dfra = 0.001d0 cc ================================================ cc > Read command line arguments: numarg = iargc ( ) do ii=1,numarg call getarg ( ii, arg(ii) ) end do cc ================================================ cc > Display help: if(arg(1).eq.'-help') then write(*,*) '' write(*,*) '_________________________________________________' write(*,*) '' write(*,*) ' LSFD' write(*,*) ' Least Squares Fourier Decomposer' write(*,*) '_________________________________________________' write(*,*) '' write(*,*) 'Usage: lsfd.e [v]' write(*,*) 'Where:' write(*,*) ' contains the time series in the form:' write(*,*) ' HJD(1) MAG(1)' write(*,*) ' HJD(1) MAG(2)' write(*,*) ' ...' write(*,*) ' HJD(N) MAG(N)' write(*,*) '

is the period in days if value is positive' write(*,*) ' is the frequency in c/d if value is negative' write(*,*) ' is the order of the Fourier fit' write(*,*) ' [v] is for verbose mode' write(*,*) '_________________________________________________' write(*,*) '' write(*,*) ' Form of the Fourier decomposition:' write(*,*) ' A_0+A_1*sin(2*pi*(t(i)-Epoch)*1/Period+Phi_1)' write(*,*) ' A_2*sin(2*pi*(t(i)-Epoch)*2/Period+Phi_2)+' write(*,*) ' A_3*sin(2*pi*(t(i)-Epoch)*3/Period+Phi_3)+..' write(*,*) '_________________________________________________' write(*,*) '' write(*,*) 'Parameters are output in the form:' write(*,*) ' period epoch a(0) Ndata sigma totamp a_i(0)' write(*,*) ' a(1) phi(1)' write(*,*) ' a(2) phi(2)' write(*,*) ' ...' write(*,*) ' a(Nord) phi(Nord)' write(*,*) '' write(*,*) 'Where: a(0) is the magnitude average' write(*,*) ' (zero-frequency constant)' write(*,*) ' a_i(0) is the intensity-averaged magnitude' write(*,*) '_________________________________________________' write(*,*) '' go to 9999 else if(numarg.lt.3) then write(*,*) '' write(*,*) 'Usage: lsfd [v]' write(*,*) 'More info: lsfd -help' write(*,*) '' go to 9999 end if cc ================================================ cc WORK BEGINS HERE cc ================================================ cc File containing the time series: filets = arg(1) cc > Read in period / frequency: read(arg(2),*) freq if (freq.lt.0.0d0) then freq = - freq else freq = 1.0d0 / freq end if cc > Read in harmonics: read(arg(3),*) norder cc > Verbose / Fixed format / Default output: isverb=0 ifixfm=0 if(numarg.gt.3 .and. arg(4).eq.'v') isverb=1 if(numarg.gt.3 .and. arg(4).eq.'f') ifixfm=1 cc > Read in time series: open(1, file = filets) nrows = 0 jj = 1 301 continue read(1,*,end = 3010) t(jj), x(jj) nrows = nrows + 1 jj = jj + 1 go to 301 3010 continue close(1) cc ndata - number of data points in the time series ndata = nrows call minmax(ndata,t,tmin,tmax) tot=tmax-tmin c avrg=0.0d0 c do ndtpnt=1,ndata c write(*,*) x(ndtpnt) c avrg = avrg + x(ndtpnt) c end do c avrg = avrg / ndata c write(*,*) avrg cc ================================================ if(ifitf.eq.0) goto 5700 freqin=freq call pfreq(norder,ndata,t,x,tot,dfra,freqin,sigm,freq,ampl) 5700 continue do ii=1,norder w(ii) = freq * ii end do call fourdec(norder,w,ndata,t,x,a0,a,b,xw,disp,t0,totamp,a0i) cc ================================================ if(ierror.eq.1) then write(*,*) '\n Enter an INTEGER for random number generator' write(*,*) '' read(*,*) iseed c iseed = 2000*iseed+1 call errors(norder,w,ndata,tot,dfra,t,x,a0,a,b,disp,t0,nsim, & iseed,amp,phs,sfreq) do 280 kk=1,norder do 281 it=1,nsim aa(it) = amp(it,kk) pp(it) = phs(it,kk) 281 continue c call stats(aa,nsim,aav,adi) c call stats(pp,nsim,pav,pdi) call omit1(asig,aa,nsim,aav,adi,nclip) call omit1(asig,pp,nsim,pav,pdi,nclip) call omit1(asig,sfreq,nsim,fav,fdi,nclip) ampavr(kk) = aav ampdsp(kk) = adi phsavr(kk) = pav phsdsp(kk) = pdi write(*,*) aav, adi, pav, pdi, fdi, nclip 280 continue end if cc ================================================ cc > Write out Fourier decomposition: cc > VERBOSE mode: write derived parameters on screen: if(isverb.eq.1) then ! VERBOSE OUTPUT write(*,*) '_________________________________________________' write(*,*) '' write(*,*) ' LSFD RESULTS' write(*,*) '_________________________________________________' write(*,*) '' write(*,801) 'NFREQ = ', norder write(*,802) 'EPOCH = ', t0 write(*,802) 'STDEV = ', disp write(*,802) 'A0 = ', a0 write(*,802) 'A0I = ', a0i write(*,802) 'A(tot)= ', totamp write(*,*) '' write(*,*) ' Frequency Amplitude Phase' write(*,*) ' -----------------------------------------' write(*,*) '' do kk=1,norder write(*,803) ' f(',kk,')=',w(kk),' a(',kk,')=', & a(kk)*1000.0d0,' p(',kk,')=',b(kk) end do write(*,*) '_________________________________________________' write(*,*) '' if(ierror.eq.1) then write(*,*) '' write(*,*) ' Errors ( NSIM = ',nsim,' ):' write(*,*) ' -----------------------------------------' write(*,*) '' do kk=1,norder write(*,803) ' f(',kk,')=',w(kk),' da(',kk,')=', & ampdsp(kk),' dp(',kk,')=',phsdsp(kk) end do write(*,*) '_________________________________________________' write(*,*) '' end if else if(ifixfm.eq.1) then ! FIXED FORMAT OUTPUT write(*,*) arg(1) write(*,175) 1.0d0/w(1), t0, a0, ndata, disp, totamp, a0i do jj=1,15 if(jj.gt.norder) then a(jj) = 0.0d0 b(jj) = 0.0d0 end if write(*,176) a(jj), b(jj) if(MODULO(jj,5).eq.0) write(*,*) '' end do else ! DEFAULT FORMAT OUTPUT write(*,175) 1.0d0/w(1), t0, a0, ndata, disp, totamp, a0i do jj=1,norder write(*,177) a(jj), b(jj) end do end if 801 format(3x,a,i10) 802 format(3x,a,f14.5) 803 format(1x,a,i2,a,f9.6,a,i2,a,f14.5,a,i2,a,f6.4) 175 format(1x,f12.8,1x,f14.6,1x,f7.4,1x,i6,1x,f6.4, & 1x,f6.4,1x,f7.4) 176 format(2x,f6.4,2x,f6.4,$) 177 format(1x,f6.4,1x,f6.4) cc > Write prewhitened time series if(iwhite.eq.1) then fileprw=trim(filets)//'1' open(6,file = fileprw) do ii=1,ndata write(6,*) t(ii), xw(ii) end do close(6) end if 9999 continue return end cc xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx cc xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx subroutine errors(norder,w,ndata,tot,dfra,t,x,a0,a,b,disp,t0,nsim, & iseed,amp,phs,sfreq) implicit real*8 (a-h,o-z) real*8 w(100) real*8 t(10000) real*8 x(10000) real*8 xwsim(10000) real*8 xsim(10000) real*8 a(100),b(100) real*8 asim(100), bsim(100) real*8 amp(1000,100), phs(1000,100) real*8 sfreq(1000) cc p7 = 2 * pi p7 = 8.0d0*datan(1.0d0) do 250 it=1,nsim do 260 ii=1,ndata s = a0 do 2001 kk=1,norder phase = w(kk) * ( t(ii) - t0 ) s = s + a(kk) * dsin( p7*phase + b(kk) ) 2001 continue xsim(ii) = s + disp*gauss(iseed) c write(20,*) t(ii), xsim(ii) 260 continue freqin=w(1) call pfreq(norder,ndata,t,xsim,tot,dfra,freqin,sigm,freq,ampl) sfreq(it)=freq call fourdec(norder,w,ndata,t,xsim,a0sim,asim,bsim, & xwsim,dispsim,t0,totamp,a0isim) do 270 kk=1,norder amp(it,kk) = 1000.0d0 * asim(kk) phs(it,kk) = bsim(kk) c write(15,*) bsim(kk) 270 continue write(*,*) it, dispsim 250 continue return end cc xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx cc xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx subroutine pfreq(m,n,t,x,tot,dfra,frmax,sigm,freq,ampl) c c This routine performs a Least-Squares search for the c best fitting frequency with a Fourier sum of order 'm' c c Input: c ~~~~~~ c m = order of the Fourier sum c n = number of the data c t(i) = time at the i-th data item c x(i) = value of the time-series at the i-th data item c dfra = accuracy of the frequency determination in terms c of the line width: dfra=df/(1/T) c frmax = starting frequency value in the iteration c c Output: c ~~~~~~~ c sigm = standard deviation of the residuals c freq = best frequency value c ampl = A_1 Fourier component at the best frequency c c--------------------------------------------------------------------- c implicit real*8 (a-h,o-z) c dimension t(10000),x(10000),xw(10000),w(100) dimension a(100),b(100) c nf = 20 df = 0.2d0/tot eps = dfra/tot freq = frmax sigm = 1.0d30 c 1 continue c if(df.lt.eps) go to 3 delf = 3.0*df fmin = freq - delf df = 2.0d0*delf/dfloat(nf-1) do 2 j=1,nf ! p0=1.0d0/(fmin+df*dfloat(j-1)) fre=fmin+df*dfloat(j-1) if(fre.lt.0.0d0) go to 2 do 4 ii=1,m w(ii) = fre * ii 4 continue call fourdec(m,w,n,t,x,a0,a,b,xw,disp,t0,totamp,a0i) ! call fourdec(p0,n,m,t,x,xw,sig,a00,aa,pp) if(disp.gt.sigm) go to 2 sigm=disp freq=fre ampl=a(1) ! print *,freq,ampl,sigm 2 continue c go to 1 c 3 continue c return end cc xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx cc xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx subroutine fourdec(norder,w,ndata, & t,x,a0,a,b,xw,disp,t0,totamp,a0i) cc This routine fits a linear least squares cc Fourier-sum to a time series. cc INPUT ARGUMENTS: cc ~~~~~~~~~~~~~~~~ cc norder - number of frequencies to be fitted cc ndata - number of data in the input time series cc ww(i) , i = 1...norder - frequencies to be fitted cc t(i) , i = 1...ndata - time cc x(i) , i = 1...ndata - time series cc OUTPUT ARGUMENTS: cc ~~~~~~~~~~~~~~~~~ cc a0 - zero frequency component cc a(i) , i = 1...norder - amplitude of the i-th component cc b(i) , i = 1...norder - phase of the j-th component cc disp - standard deviation of the fit cc t0 - epoch cc totamp - total (max-min) amplitude of light variation cc FOURIER DECOMPOSITION HAS THE FOLLOWING FORM: cc ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ cc a0 + a1 * sin( 2*pi*w1*(t-t0) + b1 ) + ... cc t(i) , x(i) , i=1...ndata - input data implicit real*8 (a-h,o-z) real*8 w(100) real*8 t(10000) real*8 te(10000) real*8 x(10000), v(10000) real*8 g(100,101) real*8 xs(10000), xw(10000) real*8 a(100),b(100) real*8 synthlc(10000) p7 = 8.0d0*datan(1.0d0) cc t0 - epoch cc Set epoch to the integer part cc of the first time of the time series: c t0 = dble(int(t(1))) c t0 = 0.0d0 t0 = t(1) do ii=1,ndata te(ii) = t(ii) - t0 end do m2 = 2 * norder m3 = m2 + 1 m4 = m3 + 1 call norm(ndata, norder, te, x, w, g) call matinv(norder, g) a0 = 0.0d+00 do j=1,m3 a0 = a0 + g(1,j) * g(j,m4) end do do 32 j=2,m3 s = 0.0d+00 do 33 k=1,m3 s = s + g(j,k) * g(k,m4) 33 continue a(j) = s 32 continue k = 0 do 37 j=2,m3,2 k = k + 1 j1 = j + 1 amp = dsqrt( a(j)*a(j) + a(j1)*a(j1) ) ph = a(j1) / amp ph = dacos(ph) if(a(j) .lt. 0.0d+00) then ph = p7 - ph end if b(k) = ph a(k) = amp 37 continue cc > CALCULATE INTENSITY AVERAGED MAGNITUDES cc > FROM THE FOURIER FIT do 5021 i=1,ndata v(i)=10.0d0**(-(x(i)-a0)/2.5d0) 5021 continue m2 = 2 * norder m3 = m2 + 1 m4 = m3 + 1 call norm(ndata,norder,te,v,w,g) call matinv(norder,g) a0i=0.0D+00 do 5022 j=1,m3 a0i=a0i + g(1,j) * g(j,m4) 5022 continue C.... A0I = INTENSITY AVERAGE FROM THE FOURIER FIT a0i=a0-2.5d0*DLOG10(a0i) cc > Calculate dispersion: disp = 0.0d0 dev = 0.0d0 do 2000 i=1,ndata s = a0 time = te(i) do 2001 k=1,norder phase = w(k) * time s = s + a(k) * dsin( p7*phase + b(k) ) 2001 continue xs(i) = s devi = x(i) - xs(i) xw(i) = devi disp = disp + devi*devi dd = xs(i) - a0 dev = dev + dd*dd 2000 continue disp = dsqrt(disp / dfloat(ndata - 2*norder - 1)) dev = dsqrt(dev / dfloat(ndata)) cc > Calculate total amplitudes: startt = t(1) stopt = startt + (1.0d0 / w(1)) * 2.0d0 nstep = 5000 tstep = (stopt - startt) / nstep do ii=0,nstep-1 s = a0 time = startt + ii * tstep do kk=1,norder phase = w(kk) * time s = s + a(kk) * dsin( p7*phase + b(kk) ) end do synthlc(ii+1) = s c write(20,*) time, s end do call minmax(nstep,synthlc,amin,amax) totamp = amax-amin return end cc xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx cc xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx SUBROUTINE MATINV(M,G) C MATRIX INVERTATION IMPLICIT REAL*8 (A-H,O-Z) C INPUT: G(M3,M3), OUTPUT: G(M3,M3) DIMENSION G(100,101) M3=2*M+1 DO 11 I=1,M3 W=G(I,I) Q=1.0D+00/W DO 12 K=1,M3 G(I,K)=Q*G(I,K) IF(K.NE.I) GO TO 12 G(I,K)=Q 12 CONTINUE DO 13 J=1,M3 IF(J.EQ.I) GO TO 13 TT=G(J,I) Q1=-TT/W DO 14 K=1,M3 G(J,K)=G(J,K)-TT*G(I,K) IF(K.NE.I) GO TO 14 G(J,K)=Q1 14 CONTINUE 13 CONTINUE 11 CONTINUE RETURN END cc xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx cc xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx SUBROUTINE NORM(N,M,T,X,F,G) C CALCULATION OF THE ELEMENTS OF THE NORMAL MATRIX AND R.H.S. IMPLICIT REAL*8 (A-H,O-Z) DIMENSION G(100,101),T(10000),X(10000),F(100) P7=8.0D0*DATAN(1.0D0) M3=2*M+1 M4=M3+1 L=0 G(1,1)=FLOAT(N) DO 1 J=2,M3,2 S=0.0D+00 C=0.0D+00 L=L+1 J1=J+1 FL=F(L) DO 2 I=1,N PH=T(I)*FL PH=P7*PH S=S+DSIN(PH) 2 C=C+DCOS(PH) G(1,J)=C 1 G(1,J1)=S L1=0 DO 3 J=2,M3,2 L1=L1+1 J1=J+1 L2=L1-1 FL1=F(L1) DO 4 K=J,M3,2 K1=K+1 L2=L2+1 FL2=F(L2) S1=0.0D+00 S2=0.0D+00 C1=0.0D+00 C2=0.0D+00 DO 5 I=1,N PH1=T(I)*FL1 PH1=P7*PH1 PH2=T(I)*FL2 PH2=P7*PH2 SPH1=DSIN(PH1) CPH1=DCOS(PH1) SPH2=DSIN(PH2) CPH2=DCOS(PH2) C1=C1+CPH1*CPH2 S1=S1+CPH2*SPH1 C2=C2+SPH2*CPH1 5 S2=S2+SPH1*SPH2 G(J,K)=C1 G(J,K1)=C2 G(J1,K)=S1 4 G(J1,K1)=S2 3 CONTINUE S=0.0D+00 DO 6 I=1,N 6 S=S+X(I) G(1,M4)=S L=0 DO 7 J=2,M3,2 S=0.0D+00 C=0.0D+00 L=L+1 J1=J+1 FL=F(L) DO 8 I=1,N PH=T(I)*FL PH=P7*PH S=S+X(I)*DSIN(PH) 8 C=C+X(I)*DCOS(PH) G(J,M4)=C 7 G(J1,M4)=S DO 9 J=1,M3 DO 10 K=J,M3 10 G(K,J)=G(J,K) 9 CONTINUE RETURN END cc xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx cc xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx FUNCTION GAUSS(ISEED) C C this routine generates Gaussian random variables of unit variance C implicit real*8 (a-h,o-z) C 1 CONTINUE P1=2.0D0*RANDOMNUMBER(ISEED)-1.0D0 P2=2.0D0*RANDOMNUMBER(ISEED)-1.0D0 RR=P1*P1+P2*P2 IF(RR.GE.1.0D0) GO TO 1 FAC=DSQRT(-2.0D0*DLOG(RR)/RR) GAUSS=P1*FAC C RETURN END cc xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx cc xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx FUNCTION RANDOMNUMBER(IX) IMPLICIT REAL*8 (A-H,O-Z) INTEGER A,P,IX,B15,B16,XHI,XALO,LEFTLO,FHI,K DATA A/16807/,B15/32768/,B16/65536/,P/2147483647/ XHI = IX/B16 XALO = (IX - XHI * B16) * A LEFTLO = XALO / B16 FHI = XHI * A + LEFTLO K = FHI / B15 IX = ((( XALO - LEFTLO * B16 ) - P ) + (FHI - K * B15 ) * B16) + K IF (IX.LT.0) IX = IX + P RANDOMNUMBER = DFLOAT(IX) * 4.65661287D-10 RETURN END cc xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx cc xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx subroutine stats(xx,nn,aver,disper) cc This routine computes the average and the dispersion of xx. implicit real*8 (a-h,o-z) real*8 xx(100) s = 0.0d0 do i=1,nn s = s + xx(i) end do aver = s / dfloat(nn) di = 0.0d0 do i=1,nn d = xx(i) - aver di=di+d*d end do disper = 0.0d0 if(nn.gt.1) then disper = dsqrt( di/dfloat(nn-1) ) end if return end subroutine omit1(asig,x,n,xave,xsig,nclip) c c--------------------------------------------------------------------- c c This routine computes average and standard deviation c of { x(i) ; i = 1,2, ..., n } by iteratively omitting c data which deviate more that asig*sigma c c Input parameters/arrays: asig, t(i), x(i), n c ~~~~~~~~~~~~~~~~~~~~~~~~ c c Output parameters: xave = average c ~~~~~~~~~~~~~~~~~~ xsig = standard deviation c nclip = number of data points clipped c c--------------------------------------------------------------------- c implicit real*8 (a-h,o-z) c dimension x(100),y(100) dimension in(100) c eps0 = 1.0d-15 nn = n c do 21 i=1,nn y(i) = x(i) 21 continue c 20 continue c call stats(y,nn,aver,disper) c c... Finding out if there is an outlier c if(disper.le.eps0) return c iout=0 devm=asig*disper do 10 i=1,nn dev=dabs(y(i)-aver) in(i)=1 if(dev.gt.devm) in(i)=0 10 continue c c... Omit outlier c k=0 do 1 i=1,nn if(in(i).eq.0) go to 1 k=k+1 y(k)=y(i) 1 continue if(k.eq.nn) go to 22 nn=k go to 20 22 continue nclip=n-nn call stats(y,nn,xave,xsig) c return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c subroutine minmax(n,x,xmin,xmax) c c This routine computes (min,max) of {x(i); i=1,...,n} c implicit real*8 (a-h,o-z) c dimension x(10000) c xmin = 1.0d30 xmax =-1.0d30 do 1 i=1,n xmin = dmin1(xmin,x(i)) xmax = dmax1(xmax,x(i)) 1 continue c return end c c