text:
        program s_pscorg
*******************************************************************************
*
* s_pscorg is a phase screen simulator which produces artificial fringe
* or speckle patterns and obtain integrated power spectrum or(/and)
* optionally specified bispectral components.
*
* The atmospheric disturbance is given by the Kolmogorov spectrum.
* Seeing condition, aperture configuration and the object ( at this
* point, only multiple point sources ) and the usage of the phase screen
* must be specified by the user in an ascii input file.
*
* Following is a sample input file. (free format, ! indicates a
* comment line.)
*
* ! Field of view: side of image in arcseconds
* Field of view
* 12
* ! seed for random number generator (large odd integer)
* Random
* 7597551
* ! outer scale, inner scale of the turbulence(m) - recommended values
* Turbulence scale
* 10.0  0.01
* ! Fried radius (r_zero) in m
* Fried
* 0.1
* ! central wavelength(m), fractional bandwidth, sampling points
* Wavelength
* 0.65e-6 0.016 3
* ! source: # of point sources
* !         intensity(arbitrary unit), position offset(x,y) in arcsecond
* Source
* 2
* 1.0  -0.7 -0.6
* 0.5   0.8 0.4
* ! # of circular apertures (including voids)
* ! configuraton: center of aperture(x,y), radius(m) and type
* ! 1 is an aperture,0 is a void. Voids should follow ordinary apertures.
* ! Example: annular mask
* Aperture
* 2
* 0.0 0.0 2.0  1
* 0.0 0.0 1.9  0
* ! Multiple subsets are taken out of one phase screen.
* ! Equally spaced sampling is specified by
* ! x_min, x_max, step size in x(m)
* ! y_min, y_max, step size in y(m)
* ! JLP/WARNING: x_min, y_min should be larger than the aperture radii !!!
* Phase_screen
* 1.25 8.75  0.14
* 1.25 8.75  2.5
* ! total number of frames to be processed
* Frames
* 100
* ------ as an option `no atmosphere' can be specified by ------
* ! obtaining beam pattern
* No_disturbance
* --------------------------------------------------------------
*
* This is a re-organized version of the pscsn.for developed on Tacos.
* The program is planned to be as users oriented as possible for the
* future development of Optical/IR interferometry.
*
* version 1.0    Dec 25, 1988 (Merry Christmas!)     Tadashi Nakajima
* Today's resolution : `I will have real holidays next year and live
*                       a human life.'
* version 2.0    Dec 30, 1988   Now sharp cut-offs of fluctuation power
* spectrum at kapmin and kapmax are introduced to investigatate the
* relative significance of the diffractive effect and refractive effect.
* For instance, the refractive effect can be seen by setting lo = 0.1,
* while the diffractive effect can be investigated by setting L0 = 0.4.
*
* JLP:
*   - I introduce fourn1 and jlp_random
*   - I remove the common blocks, split the program into subroutines
*      and complete the documentation
*   - AND I CORRECT DEEPLY THIS PROGRAM SINCE IT RELIED ON AN UNKNOWN
*      STRUCTURE OF FOURIER COMPLEX ARRAYS AND FOURIER TRANSFORM!
*   - I add r_zero and field of view into the parameter file
*   - I change the screen generator (jlp_phase_screen) so that the
*   simulation is compatible with the r_zero value (the old version was not)
*   - I replace the pseudo-convolution of the objet with the PSF by a true
*   convolution
*   - I implement a good formula for the image smearing due to finite bandwidth
*   (the original was wrong!)
*
* JLP
* Version 26/05/2008
*******************************************************************************
*------ declarations ---------*
        implicit none
* Constants
        real*4 pi

* Integrated structure constant of the phase fluctuation
* LCn2 = \int Cn2(h) dh
* And also:
* r_0 = (0.423 kk^2 LCn2)^{-3/5}
        real*4 LCn2, r_zero

* Size of phase screen and resolution (m/pix)
        real*4 Ls, Lp

* Size of phase screen, size of power spectrum
        integer*4 Ns, Np
        parameter(Ns = 512)
        parameter(Np = Ns/2)

* Wave number, scale length of turbulence
* fov: field of view (for the simulated images)
        real kk, L0, lambda, lo, fov

* psc: complex phase screen, real and imag are independent.
* FT_object(Np,Np): Fourier Transform of the object
        complex psc(Ns,Ns), FT_object(Np,Np)

* Power spectrum array
        real*4 pwr(Np,Np)

* Input control parameter file and output power spectrum file
        character ifile*60, ofile*60, ocomments*80

* I/O units
        integer*4 f_out
        parameter(f_out = 12)

* nframes: number of frames
* nscreens: number of phase screens
* frm_per_scr: number of frames per phase screen
        integer*4  nframes, nscreens, frm_per_scr

* Maxcomp: maximum number of components in source
* Maxap: maximum number of apertures
        integer*4 Maxcomp, Maxap
        parameter(Maxcomp = 50, Maxap = 15)

* Source information (intensity and position in arcsec)
        real*4 s_int(Maxcomp), s_pos(Maxcomp,2), object(Np,Np)

* Aperture information (coordinate, radius, type)
        real*4 c_pos(Maxap,2), c_rad(Maxap)
        integer*4 c_typ(Maxap)
* Miscellaneous:
        real*4 zero, x_min, x_max, xstep, y_min, y_max, ystep
        real*4 rapsec, fractional_bw
        integer*4 ncalls, ncircle, nsamp, lx, ly, i, iz, jz, iframe
        integer*4 ncomp, seed, no_dist

* Initialization of variables:
        pi = 3.141592653

* rapsec: radians per arcsecond
        rapsec = pi/(180. * 3600.0)

* JLP
* Program can then be run as "runs s_pscorg s_screen.input power"
        call jlp_begin
        call jlp_inquifmt

* Input file names (from terminal or "jlp_lu5.tmp" file when run with "runs")
        write(6,10)
   10   format(1x,'Input ascii file name : ',$)
        read(5,'(a)') ifile
        write(6,20)
   20   format(1x,'Output power FITS file name : ',$)
        read(5,'(a)') ofile
        write(6,25)
   25   format(1x)

C Read input parameter file:
        call read_input_file(ifile, seed, L0, lo, lambda, fractional_bw,
     1                nsamp, ncomp, s_int, s_pos, ncircle, c_pos,
     1                c_rad, c_typ, x_min, x_max, xstep, y_min,
     1                y_max, ystep, fov, r_zero, nframes, no_dist,
     1                Maxcomp, Maxap)

C Initialization of random generator:
        call jlp_random_init(seed)

*  kk: wave number
        kk = 2.0 * pi / lambda

* Original version with LCn2 = 5.0e-13
* r_0 = (0.423 kk^2 LCn2)^{-3/5}
        LCn2 = (r_zero**(3./5.)) / (0.423 * kk * kk)
        write(6,287) lambda, r_zero, LCn2
287     format('lambda=',e12.5,' r_zero=',f8.3,' LCn2=',e12.5)

* Scale of image in pixels/arcsecond: Scale = (Ls/lambda)*(Np/Ns)*rapsec
* Field of view = Np/Scale
* Size of phase screen and resolution (m/pix)
        Ls = Ns / ((fov / lambda) * rapsec)
        write(6,*) 'Size of pupil: Ls = ',Ls,' m'
        Lp = Ls / real(Ns)
        write(6,*) 'Resolution of pupil: Lp = ',Lp,' m'

C Create an image of the source from ncomp, s_int, s_pos
         call create_object(object, FT_object, Np, Ns, ncomp,
     1              s_int, s_pos, Ls, lambda, Maxcomp)

* Number of frames per phase screen
        frm_per_scr =
     1  (int((x_max-x_min)/xstep)+1)*(int((y_max-y_min)/ystep)+1)
*
        if( mod(nframes,frm_per_scr) .eq. 0 ) then
          nscreens = nframes / frm_per_scr
        else
          nscreens = nframes / frm_per_scr + 1
        end if

        do iz=1,Np
        do jz=1,Np
          pwr(jz,jz) = 0.
        end do
        end do
* Two real phase screens per complex random number generation in Fourier
* domain.

        if( mod(nscreens,2) .eq. 0 ) then
          ncalls = nscreens/2
        else
          ncalls = nscreens/2 + 1
        end if
        write(6,*)' ******** Will now simulate ', nframes,
     1            ' elementary frames ************'
        write(6,*)' Number of frames per phase screen: ', frm_per_scr
        write(6,*)' Number of phase screen needed: ', nscreens
        write(6,*)' Number of calls of phase_screen routine needed: ',
     1            ncalls

        iframe = 0
        do 300 i = 1, ncalls
*------------------------------------
        if (no_dist.eq.1) then
* no atmospheric disturbance
        write(6,312) i, iframe+1
  312   format('+', 'generating ',i3,
     1         'th zero phase screen (iframe=',i3,')')
          do iz=1, Ns
            do jz=1, Ns
              psc(iz,jz) = (0.0,0.0)
            end do
          end do
        else
* Create a complex phase screen
        write(6,310) i, iframe+1
  310   format('+', 'generating ',i3,
     1         'th complex phase screen (iframe=',i3,')')
C New version:
        call jlp_phase_screen(psc, Ns, seed, L0, r_zero)
C Old version:
C        call  phase_screen(psc, Ns, seed, Ls, L0, lo,
C     1                     kk, LCn2)
        end if

* sample phase screen, form images, calculate power
        call  power(psc, Ns, FT_object, pwr, Np, iframe, nframes, Ls,
     1                  lambda, fractional_bw, nsamp, ncomp,
     1                  s_int, s_pos, ncircle, c_pos, c_rad, c_typ,
     1                  x_min, x_max, xstep, y_min, y_max, ystep,
     1                  Maxcomp, Maxap)

          if(iframe .ge. nframes) goto 400
*
*------------------------------------
* End of loop
  300   continue
*

* Exit from loop
  400   continue

        call jlp_recent_real(pwr, Np)
* Normalize integrated power spectrum
        zero = pwr(Np/2+1,Np/2+1)
        do lx=1, Np
          do ly=1, Np
            pwr(lx,ly) = pwr(lx,ly)/zero
          end do
        end do

* Create a FITS file
        write(6,*)' Writing output file:',ofile
        write(ocomments,88)ifile,char(0)
88      format('From ',A,A)
        call jlp_writeimag(pwr,Np,Np,Np,ofile,ocomments)
*
        call jlp_end
        stop
        end
***************************************************************************
* Creates complex phase screen (JLP's version)
*
* INPUT:
* Ns: size of phase screen (512)
* seed: seed for the random number generator,
*
* OUPUT:
*
***************************************************************************
        subroutine jlp_phase_screen(psc, Ns, seed, L0, r_zero)
        implicit none
        integer*4 Ns, seed
        complex psc(Ns,Ns),iiu
        real*4 L0, r_zero, pi, work(Ns,Ns), mean, sigma
        real*4 rad2
        real*8 sum, sumsq
        logical debug
        character out_file*60, out_comments*80
        real*4 urad, radi, Fs, Fs0, utheta

* N(): dimension array for fft
* Ns_cent: origin of the phase screen
        integer*4 N(2), Ns_cent, jx, jy, i, j

        debug=.true.
* Constants:
* iiu = sqrt(-1) or "i"
        iiu = (0.0,1.0)
        pi = 3.141592653

* r_0 = (0.423 kk^2 LCn2)^{-3/5}
* Hence Fs0 is proportional to kk^2
C JLP2008 / I put Ns^4 instead of Ns^2 ...
        Fs0 = 0.023 * Ns * Ns * Ns * Ns * (L0 / r_zero)**(5./3.)

*------ generation of random phase in spatial domain -------------
* Complex Gaussian random numbers are generated by
* the inverse transform method.

        N(1) = Ns
        N(2) = Ns
        Ns_cent = Ns/2+1

        do 100 jx=1, Ns
          do 110 jy=1, Ns

            rad2=(jx-Ns_cent)*(jx-Ns_cent)+(jy-Ns_cent)*(Jy-Ns_cent)

* rad**(-11/3) is undefined at rad2=0.
            if(rad2 .eq. 0.0) rad2 = 0.5

10          call jlp_random(urad)
            if ( urad.gt.0.0 ) then
              radi = sqrt(-2.0*log(urad))
            else
              goto 10
            end if
            Fs = Fs0 * rad2**(-11.0/6.0)
            call jlp_random(utheta)
* iiu = sqrt(-1) or "i"
* Hence psc is proportional to kk = 2 pi / lambda
            psc(jx,jy) = sqrt(Fs)*radi*cexp(iiu*2.0*pi*utheta)
  110     continue
  100   continue

*---------- Fourier transform ------------------------------
* Shift origin of Fourier space to (0,0) before inverse Fourier transform
* Old recentering version from original program was wrong!
        call jlp_recent_cmplx(psc, Ns)

        write(6,555) N(1),N(2)
  555   format('+', 'start ',i3,'x',i3,' sized fft              ')
C        call fourn1(psc,N,2,-1,+1)
C JLP2008 (by comparing with s_psc.for):
C fourn1(data,nn,ndim,isign)
C SHOULD USE MY fourn1 routine (in fft_jlp.for)
C NOT s_fourn1.for !!!
        call fourn1(psc,N,2,-1)
        write(6,666)
  666   format('+', 'end fft for screen                         ')

* Now psc is a complex random phase screen
* to avoid discontinuity shift origin of phase screen to (Ns/2,Ns/2)
* Old recentering version from original program was wrong!
        call jlp_recent_cmplx(psc, Ns)

C Output phase screen:
        if(debug)then
          sum = 0.
          sumsq = 0.
          do i=1, Ns
            do j=1, Ns
              work(i, j) = real(psc(i,j))
              sum = sum + work(i,j)
              sumsq = sumsq + work(i,j) * work(i,j)
            end do
          end do
          mean = sum / real(Ns * Ns)
          sigma = sqrt(sumsq / real(Ns * Ns) - mean * mean)
          write(6,*)'Phase screen: mean=', mean,' sigma=',sigma
          out_file='screen_real_psc'
          out_comments='Real part of phase screen'
          call jlp_writeimag(work,Ns,Ns,Ns,out_file,out_comments)
        endif

        return
        end
***************************************************************************
* Creates complex phase screen (original version)
*
* INPUT:
* Ns: size of phase screen (512)
* seed: seed for the random number generator,
* kk: wave number
* LCn2 integrated CN2 over all the altitudes (originally: LCn2=5.0e-13)
*
* OUPUT:
*
***************************************************************************
        subroutine phase_screen(psc, Ns, seed, Ls, L0, lo,
     1                          kk, LCn2)
        implicit none
        integer*4 Ns, seed
        complex psc(Ns,Ns),iiu
        real*4 dkap, kapmin, kapmax, kk, LCn2, pi, Ls, L0, lo, Lp
        real*4 work(Ns,Ns), mean, sigma, ww
        real*8 sum,sumsq
        logical debug
        character out_file*60, out_comments*80

* kap: radial spatial frequency in aperture plane ( * 2 pi / L)
        real*4 kap, atten, urad, radi, Fs, Fs0, utheta

* N(): dimension array for fft
* Ns_cent: origin of the phase screen
        integer*4 N(2), Ns_cent, jx, jy, i, j

        debug=.true.
* Constants:
* iiu = sqrt(-1) or "i"
        iiu = (0.0,1.0)
        pi = 3.141592653

* Lp: pixel scale of phase screen
        Lp = Ls/Ns

* kapmin, kapmax: min and max spatial frequencies in aperture plane
* Ls: size of phase screen (in m)
* Lp: resolution, i.e., size of one pixel in pupil (in m/pix)
        if (Ls .lt. L0) then
          kapmin = 2.0 * pi / Ls
        else
          kapmin = 2.0 * pi / L0
        end if

        if (lo .lt. Lp) then
           kapmax = 2.0 * pi / Lp
        else
           kapmax = 2.0 * pi / lo
        end if

* dkap: one pixel in spatial frequency in aperture plane (dkap = 2 pi / Ls)
        dkap = 2.0 * pi / Ls

* dkap: one pixel in spatial frequency in aperture plane (dkap = 2 pi / Ls)
* r_0 = (0.423 kk^2 LCn2)^{-3/5}
        Fs0 = 0.033*pi*LCn2*kk*kk*2.0*dkap*dkap

*------ generation of random phase in spatial domain -------------
* Complex Gaussian random numbers are generated by
* the inverse transform method.

        N(1) = Ns
        N(2) = Ns
        Ns_cent = Ns/2+1

        do 100 jx=1, Ns
          do 110 jy=1, Ns

            kap = dkap*sqrt( float( (jx-Ns_cent)**2+(jy-Ns_cent)**2 ) )

* Introducing smooth cut-offs below kapmin and beyond kapmax.

            if(kap .gt. kapmax) atten = exp(-(kap-kapmax)/(5.0*dkap))
            if(kap .lt. kapmin) atten = exp(-(kapmin-kap)/(5.0*dkap))
            if(kap.le.kapmax .and. kap.ge.kapmin) atten = 1.0

* kap**(-11/3) is undefined at kap=0.
            if(kap .eq. 0.0) kap = dkap

10          call jlp_random(urad)
            if ( urad.gt.0.0 ) then
              radi = sqrt(-2.0*log(urad))
            else
              goto 10
            end if
* Fs0 = 0.033*pi*LCn2*k*k*2.0*dkap*dkap
* Fs = 0.033*pi*LCn2*kk*kk*2.0*kap**(-11.0/3.0)*dkap*dkap
            Fs = Fs0 * kap**(-11.0/3.0)
            call jlp_random(utheta)
* iiu = sqrt(-1) or "i"
            psc(jx,jy) = atten*sqrt(Fs)*radi*cexp(iiu*2.0*pi*utheta)
  110     continue
  100   continue

*---------- Fourier transform ------------------------------
* Shift origin of Fourier space to (0,0)
* Old recentering version from original program was wrong!
        call jlp_recent_cmplx(psc, Ns)

        write(6,555) N(1),N(2)
  555   format('+', 'start ',i3,'x',i3,' sized fft              ')
C        call fourn1(psc,N,2,-1,+1)
C JLP2008 (by comparing with s_psc.for):
C fourn1(data,nn,ndim,isign)
C SHOULD USE MY fourn1 routine (in fft_jlp.for)
C NOT s_fourn1.for !!!
        call fourn1(psc,N,2,-1)
        write(6,666)
  666   format('+', 'end fft for screen                         ')

C JLP2008: I multiply with Ns*Ns, since otherwise
C the values of psc are much too small!
        ww = real(Ns * Ns)
          do i=1, Ns
            do j=1, Ns
              psc(i,j) = psc(i,j) * ww
            end do
          end do

* Now psc is a complex random phase screen
* to avoid discontinuity shift origin of phase screen to (Ns/2,Ns/2)
* Old recentering version from original program was wrong!
        call jlp_recent_cmplx(psc, Ns)

C Output phase screen:
        if(debug)then
          sum = 0.
          sumsq = 0.
          do i=1, Ns
            do j=1, Ns
              work(i, j) = real(psc(i,j))
              sum = sum + work(i,j)
              sumsq = sumsq + work(i,j) * work(i,j)
            end do
          end do
          mean = sum / real(Ns * Ns)
          sigma = sqrt(sumsq / real(Ns * Ns) - mean * mean)
          write(6,*)'Phase screen: mean=', mean,' sigma=',sigma
          out_file='screen_real_psc'
          out_comments='Real part of phase screen'
          call jlp_writeimag(work,Ns,Ns,Ns,out_file,out_comments)
        endif

        return
        end
*******************************************************************************
* power
* Produces fringe/speckle patterns of multiple point sources.
* Finite bandwidth effect is taken into account.
*
* INPUT:
* psc(Ns,Ns): phase screen
* Ns: size of phase screen (=512)
* Np: size of power spectrum (=Ns/2)
* iframe: index of current frame to be processed
* nframes: number of frames to be simulated
* ncomp: number of components for the source
* s_int, s_pos: source information (intensity and position in arcsec)
* ncirc: number of circles for the aperture
* c_pos, c_rad, c_typ: aperture information (coordinate, radius, type)
* Ls: size of phase screen (in m)
* lambda, fractional_bw, nsamp: Wavelength parameters: central wavel.,
*                          fractional bandwidth and sampling points
* x_min, x_max, xstep : sampling parameters to extract a phase screen
*                       from a larger one
* y_min, y_max, ystep : sampling parameters to extract a phase screen
*                       from a larger one
* Maxcomp: maximum number of components in source
* Maxap: maximum number of apertures
*
* OUTPUT:
* pwr(Np,Np): power spectrum
*******************************************************************************
        subroutine power(psc, Ns, FT_object, pwr, Np, iframe,
     1                  nframes, Ls,
     1                  lambda, fractional_bw, nsamp, ncomp,
     1                  s_int, s_pos, ncircle, c_pos, c_rad, c_typ,
     1                  x_min, x_max, xstep, y_min, y_max, ystep,
     1                  Maxcomp, Maxap)
        implicit none
        integer Ns, Np, Maxcomp, Maxap
        complex psc(Ns,Ns), FT_object(Np,Np)
        real*4 pwr(Np,Np), Ls, lambda, fractional_bw
        real*4 x_min, x_max, xstep, y_min, y_max, ystep
        integer*4 iframe, nframes, nsamp, ncomp, ncircle

* Aperture information (coordinate, radius, type)
        real*4 c_pos(Maxap,2), c_rad(Maxap)
        integer*4 c_typ(Maxap)

* Source information (intensity and position in arcsec)
        real*4 s_int(Maxcomp), s_pos(Maxcomp,2)

* Temporary arrays
        real*4  PSF1(Np,Np), PSF2(Np,Np)

C Miscellaneous:
        integer*4 N(2), Np_cent
        integer*4 output_screen
        real*4  Lp, Scale, x, y

* Constants
        real*4 pi, rapdeg, rapsec
        complex iiu
* iiu = sqrt(-1) or "i"
        iiu = (0.0,1.0)
        pi = 3.141592653

* rapsec: radians per arcsecond
        rapdeg = pi/180.0
        rapsec = rapdeg/3600.0

        output_screen = 0
        N(1) = Np
        N(2) = Np
* Ls: size of phase screen (in m)
* Lp: pixel scale of phase screen
        Lp = Ls/Ns
        Np_cent = Np/2+1
* lambda: wavelength (in meters)
* Scale: in pixels/arcsecond
        Scale = (Ls/lambda)*(real(Np)/real(Ns))*rapsec

* Setting apertures on Np sized window in Ns sized complex phase screen

        do 110 x=x_min,x_max,xstep
        do 120 y=y_min,y_max,ystep
**************************** major loop ***************************************

        write(6,111) x,y
  111   format('+',' creating a pair of speckle/fringe patterns at',
     1         'x=',f10.3,' y=',f10.3)

        write(6,333) ncircle, nsamp, c_typ(1), c_typ(2)
  333   format(1x,'ncircle nsamp typ(1) typ(2) ',1x,4(i2,1x) )

C Compute two elementary PSFs (PSF1 and PSF2)
        call compute_elementary_PSFs(psc, x, y, Ls,
     1                  lambda, fractional_bw, nsamp,
     1                  ncircle, c_pos, c_rad, c_typ,
     1                  PSF1, PSF2, Ns, Np, Maxap, output_screen)

C Convolution of the object with the elementary PSFs

        write(6,320) iframe+1
  320   format('+', 'processing ',i5,' th frame                    ')
       call process_elementary_image(pwr, PSF1, FT_object,
     1                      Np, output_screen)
        iframe = iframe + 1
        if( iframe .ge. nframes) return

        write(6,320) iframe+1
       call process_elementary_image(pwr, PSF2, FT_object,
     1                      Np, output_screen)
        iframe = iframe + 1
        if( iframe .ge. nframes) return

************************* end of major loop ***********************************
  120   continue
  110   continue

        return
        end
***************************************************************************
* read_input_file: read input parameter file
*
* INPUT:
* ifile: name of input parameter file
*
* OUTPUT:
* seed: seed for random number generator (large odd integer)
* L0, lo: external and internal scale lengths of turbulence
* lambda, fractional_bw, nsamp: Wavelength parameters: central wavel.,
*                          fractional bandwidth and sampling points
* x_min, x_max, xstep : sampling parameters to extract a phase screen
*                       from a larger one
* y_min, y_max, ystep : sampling parameters to extract a phase screen
*                       from a larger one
* ncomp: number of components for the source
* s_int, s_pos: source information (intensity and position in arcsec)
* ncirc: number of circles for the aperture
* c_pos, c_rad, c_typ: aperture information (coordinate, radius, type)
* fov: field of view (in arcseconds)
* nframes: number of frames to be simulated
* no_dist: flag set to one if no atmospheric disturbance
* Maxcomp: maximum number of components in source
* Maxap: maximum number of apertures
*
***************************************************************************
        subroutine read_input_file(ifile, seed, L0, lo,
     1                  lambda, fractional_bw, nsamp, ncomp,
     1                  s_int, s_pos, ncircle, c_pos, c_rad, c_typ,
     1                  x_min, x_max, xstep, y_min, y_max, ystep,
     1                  fov, r_zero, nframes, no_dist, Maxcomp, Maxap)
        implicit none
        integer*4 Maxcomp, Maxap
        character ifile*60, comment*80
        real*4 L0, lo, lambda, fractional_bw
        real*4 s_int(Maxcomp), s_pos(Maxcomp,2)
        real*4 c_pos(Maxap,2), c_rad(Maxap)
        real*4 x_min, x_max, xstep, y_min, y_max, ystep, fov, r_zero
        integer*4 seed, nsamp, ncomp, ncircle, nframes, no_dist
        integer*4 c_typ(Maxap), j, n
* I/O units
        integer*4 f_in
        parameter(f_in = 11)

* Initialize no_dist to default value:
        no_dist = 0

* Read the input file
        open(unit=f_in, file=ifile, status='old')

   30   format(a80)
        do 40 j=1,100

          read(f_in,30,end=100) comment
          if( comment(1:1) .eq. '!') goto 40

          write(6,30) comment

C Keyword "Field of view", parameter: fov
          if(comment(1:5) .eq. 'Field')then
            read(f_in,*) fov
            write(6,*) 'OK: fov=',fov,' arcseconds'
C Keyword "Random", parameter: seed (random generator):
          else if(comment(1:6) .eq. 'Random')then
            read(f_in,*) seed
            write(6,*) 'OK: seed=',seed
C Keyword "Turbulence scale", parameters: outer and inner  scales
          else if(comment(1:10) .eq. 'Turbulence')then
            read(f_in,*) L0, lo
            write(6,*) 'OK: L0=',L0,' lo=',lo
C Keyword "Fried radius", parameters: outer and inner  scales
          else if(comment(1:5) .eq. 'Fried')then
            read(f_in,*) r_zero
            write(6,*) 'OK: r0=',r_zero
C Keyword "Wavelength", parameters: central wavel., fractional bandwidth,
C sampling points
C lambda: wavelength (in meters)
C Fractional bandwidth = 2 (Fh - Fl)/(Fh + Fl)
C where Fh is the highest frequency limit with signal 10dB below peak emission
C where Fl is the lowest frequency limit with signal 10dB below peak emission
          else if(comment(1:10) .eq. 'Wavelength')then
            read(f_in,*) lambda, fractional_bw, nsamp
            write(6,*) 'OK: lambda=',lambda,
     1           'm  fractional_bw=',fractional_bw,' nsamp=', nsamp,
     1           '(hence bandwidth =',fractional_bw * lambda,' m)'
C Keyword "Source", parameters: number of points, followed by
C Intensity1, X1, Y1
C Intensity2, X2, Y2...
          else if(comment(1:6) .eq. 'Source')then
            read(f_in,*) ncomp
            write(6,*) 'OK: Source with ',ncomp,' components'
            do n=1,ncomp
              read(f_in,*) s_int(n),s_pos(n,1),s_pos(n,2)
              write(6,*)' int, x, y: ',s_int(n),s_pos(n,1),s_pos(n,2)
            end do
C Keyword "Aperture", parameters: number of apertures, followed by
C Xcenter, Ycenter, radius, type (0 or 1)
          else if(comment(1:8) .eq. 'Aperture')then
            read(f_in,*) ncircle
            write(6,*) 'OK: Aperture with ',ncircle,' circles'
            do n=1,ncircle
              read(f_in,*) c_pos(n,1),c_pos(n,2),c_rad(n),c_typ(n)
              write(6,*)' x, y, rad, type: ',
     1                    c_pos(n,1),c_pos(n,2),c_rad(n),c_typ(n)
            end do
C Keyword "Phase screen", parameters: x_min, x_max, step
          else if(comment(1:5) .eq. 'Phase')then
            read(f_in,*) x_min, x_max, xstep
            read(f_in,*) y_min, y_max, ystep
            write(6,*) 'OK: xmin, xmax, xstep:', x_min, x_max, xstep
            write(6,*) 'OK: ymin, ymax, ystep:', y_min, y_max, ystep
C Keyword "Frame", parameters: nber of frames
          else if(comment(1:6) .eq. 'Frames')then
            read(f_in,*) nframes
            write(6,*) 'OK: nframes=',nframes
C Keyword "No disturbance" (no atmosphere)
          else if( comment(1:2) .eq. 'No') then
            no_dist = 1
            write(6,*) 'OK: no atmospheric disturbance'
          end if
C EOF loop
   40   continue

C Switch to label=100 in case of end of file:
  100   close(unit=f_in)
        return
        end
*****************************************************************
* Recenter Fourier transform
* and shift origin from (0,0) to (Ns/2+1, Ns/2+1) or inversely
*
* Old recentering version from original program was not compatible
* with "fourn1" (and relative to an unknown data structure)
*
* Quadrants
*     3 4        2 1
*     1 2    to  4 3
*****************************************************************
       subroutine jlp_recent_cmplx(psc, Ns)
       implicit none
       integer Ns
       complex psc(Ns,Ns), tmp
       integer i, j
C Quadrants
C     3 4        2 1
C     1 2    to  4 3
        do i=1,Ns/2
          do j=1,Ns/2
C Flip quadrants 2 and 3:
            tmp = psc(i+Ns/2,j)
            psc(i+Ns/2,j) = psc(i,j+Ns/2)
            psc(i,j+Ns/2) = tmp
C Flip quadrants 1 and 4:
            tmp = psc(i+Ns/2,j+Ns/2)
            psc(i+Ns/2,j+Ns/2) = psc(i,j)
            psc(i,j) = tmp
          end do
        end do

        return
        end
*****************************************************************
* Recenter Fourier transform
* and shift origin from (0,0) to (Ns/2+1, Ns/2+1) or inversely
*
* Old recentering version from original program was not compatible
* with "fourn1" (and relative to an unknown data structure)
*
* Quadrants
*     3 4        2 1
*     1 2    to  4 3
*****************************************************************
       subroutine jlp_recent_real(pwr, Np)
       implicit none
       integer Np
       real pwr(Np,Np), tmp
       integer i, j
C Quadrants
C     3 4        2 1
C     1 2    to  4 3
        do i=1,Np/2
          do j=1,Np/2
C Flip quadrants 2 and 3:
            tmp = pwr(i+Np/2,j)
            pwr(i+Np/2,j) = pwr(i,j+Np/2)
            pwr(i,j+Np/2) = tmp
C Flip quadrants 1 and 4:
            tmp = pwr(i+Np/2,j+Np/2)
            pwr(i+Np/2,j+Np/2) = pwr(i,j)
            pwr(i,j) = tmp
          end do
        end do

        return
        end
***********************************************************************
* To create an image corresponding to the source
*
* INPUT:
* Np: size of object array (=Ns/2)
* Ns: size of phase screen (512)
* ncomp: number of components for the source
* s_int, s_pos: source information (intensity and position in arcsec)
* Ls: size of phase screen (in m)
* lambda: wavelength (in meters)
* Maxcomp: maximum number of components in source
*
* OUTPUT:
* object(Np,Np): array with ones at the location of the components
* FT_object(Np,Np): Fourier Transform of the object
***********************************************************************
        subroutine create_object(object, FT_object, Np, Ns, ncomp,
     1                          s_int, s_pos, Ls, lambda, Maxcomp)
        implicit none
        integer Np, Ns, ncomp, Maxcomp
        real*4 s_int(Maxcomp), s_pos(Maxcomp,2)
        real*4 object(Np,Np), Ls, lambda, Scale, pi, rapsec
        complex FT_object(Np,Np)
        integer i, j, ix, iy, Ncent, N(2)
        character out_file*60, out_comments*80

* Initialization of variables:
        pi = 3.141592653

* rapsec: radians per arcsecond
        rapsec = pi/(180. * 3600.0)

C Initializing the array:
        do i=1,Np
          do j=1,Np
            object(i,j) = 0.
          end do
        end do

* Scale: in pixels/arcsecond
        Scale = (Ls/lambda)*(real(Np)/real(Ns))*rapsec
        Ncent = Np/2 + 1
        write(6,*)' Scale:', Scale,' pixels/arcsec'

C Loop on all components:
        do i=1,ncomp
* ncomp components in source structure
          ix = Ncent + int(Scale * s_pos(i,1))
          iy = Ncent + int(Scale * s_pos(i,2))
          ix = min(Np, max(1, ix))
          iy = min(Np, max(1, iy))
          write(6,*)' Component #',i,' ix=',ix,' iy=',iy
          object(ix,iy) = object(ix,iy) + s_int(i)
        end do

C Transfer object to FT_object for Fourier transform
        do i=1,Np
          do j=1,Np
            FT_object(i,j)  = complex(object(i,j), 0.)
          end do
        end do

        N(1) = Np
        N(2) = Np
        call fourn1(FT_object,N,2,1)

C Output object for a diagnostic:
          out_file='s_object'
          out_comments='Source'
          call jlp_writeimag(object,Np,Np,Np,out_file,out_comments)
        return
        end
**************************************************************************
* Compute two elementary PSFs PSF1 and PSF2
* from a subwindow of the phase screen starting at x,y
* Take into account wavelength smearing due to finite bandpass
*
* INPUT:
* psc: phase screen to be used for computing the elementary PSFs
* lambda: wavelength (in meters)
* fractional_bw: fractional bandwidth
* nsamp: number of sampling points to take the finite bandwidth into account
* x, y : origin of the subwindow used to extract a phase screen from psc
* ncirc: number of circles for the aperture
* c_pos, c_rad, c_typ: aperture information (coordinate, radius, type)
*
* OUTPUT:
* PSF1, PSF2: elementary PSFs corresponding to real and imaginary
*                  parts of the phase screen
**************************************************************************
        subroutine compute_elementary_PSFs(psc, x, y, Ls,
     1                  lambda, fractional_bw, nsamp,
     1                  ncircle, c_pos, c_rad, c_typ,
     1                  PSF1, PSF2, Ns, Np, Maxap, output_screen)
        implicit none
        real*4 x, y, Ls, lambda, fractional_bw
        integer Ns, Np, Maxap
        complex psc(Ns,Ns)
        real*4  PSF1(Np,Np), PSF2(Np,Np)
        integer nsamp, ncircle, output_screen
* Aperture information (coordinate, radius, type)
        real*4 c_pos(Maxap,2), c_rad(Maxap)
        integer*4 c_typ(Maxap)

* Temporary arrays
        complex PupilFunc1(Np,Np), PupilFunc2(Np,Np)

* Miscellaneous:
        real work(Np,Np), ww
        real factor, xa, ya, ox, oy, Lp
        integer loop1, jx, jy, ix, iy, lx, ly, i, Np_cent
        integer iz, jz, N(2)
        real Delta_lambda, lambda_min, l_step
        character out_file*60, out_comments*80
        complex iiu
* iiu = sqrt(-1) or "i"
        iiu = (0.0,1.0)

        Lp = Ls/Ns
        Np_cent = Np/2+1

C Reset image arrays
        do jx=1,Np
          do jy=1,Np
            PSF1(jx,jy) = 0.0
            PSF2(jx,jy) = 0.0
          end do
        end do

* Maximum variation of lambda is Delta_lambda
        Delta_lambda = fractional_bw * lambda
        lambda_min = lambda * ( 1. - fractional_bw / 2.)
        l_step = Delta_lambda / nsamp

******************* BEGIN OF LOOP1 ********************
********* on all the (nsamp) wavelengths belonging to the bandwidth *******
        do 100 loop1=1,nsamp
**************** intermediate loop *****************
* nsamp data points within the bandpass
* Old version: wrong!
*         factor =  1.0 + fractional_bw*(loop1 - INT(nsamp/2))
* JLP's version since psc is proportional to kk = 2 pi / lambda:
         factor =  lambda / (lambda_min + loop1 * l_step)
         write(6,*)'loop1=',loop1,' factor=',factor

C Reset temporary arrays
        do jx=1,Np
          do jy=1,Np
             PupilFunc1(jx,jy) = (0.0,0.0)
             PupilFunc2(jx,jy) = (0.0,0.0)
          end do
        end do

C Loop on all the circular apertures
          do 30 i=1,ncircle
C FOR DEBUG:
C            write(6,*)'Circle #',i,' centered on x+cposx, y+cposy'
C            write(6,*) '  x=',x,' y=',y,' cposx=',c_pos(i,1)
C            write(6,*)' c_posy=',c_pos(i,2),' c_rad(i)=',c_rad(i)
C
*-----------------------------------------
*       filling pdata with aperture phases

           do 40 xa = x+c_pos(i,1)-c_rad(i), x+c_pos(i,1)+c_rad(i), Lp
           do 50 ya = y+c_pos(i,2)-c_rad(i), y+c_pos(i,2)+c_rad(i), Lp
             ox = xa-x-c_pos(i,1)
             oy = ya-y-c_pos(i,2)
             if( ox*ox+oy*oy .lt. c_rad(i)*c_rad(i) ) then
* Np_cent = Np/2+1
               ix = nint( (xa-x)/Lp )+Np_cent
               iy = nint( (ya-y)/Lp )+Np_cent

* Shift origin of the window to (Np_cent,Np_cent)
* Lp: pixel scale of phase screen (Lp = Ls/Ns)
               lx = nint(xa/Lp)
               ly = nint(ya/Lp)
               if(lx.lt.0.or.ly.lt.0) then
                write(6,*)'Fatal error: lx=',lx,' ly=',ly
                write(6,*)'Check if xmin > radius in parameter file...'
                stop
               endif

* iiu = sqrt(-1) or "i"
* NB: psc is proportional to kk = 2 pi / lambda
* correction by factor = lambda / (lambda + dlambda):
                PupilFunc1(ix,iy) = cexp(iiu*factor
     1                                *real(psc(lx,ly)) )*c_typ(i)
                PupilFunc2(ix,iy) = cexp(iiu*factor
     1                                *aimag(psc(lx,ly)) )*c_typ(i)

              end if
   50       continue
   40       continue
*-----------------------------------------
   30     continue

C Output real part of complex pupil function:
        if(output_screen.eq.0)then
          do iz=1, Np
            do jz=1, Np
              work(iz, jz) = real(PupilFunc1(iz,jz))
            end do
          end do
          out_file='s_PupilFunc1'
          out_comments='real part of exp (i phase screen)'
          call jlp_writeimag(work,Np,Np,Np,out_file,out_comments)
          output_screen = 1
        endif

* Simulate propagation from aperture plane to detector plane

C        call fourn1(PupilFunc1,N,2,-1,+1)
C        call fourn1(PupilFunc2,N,2,-1,+1)
C JLP2008 (by comparing with s_psc.for):
C fourn1(data,nn,ndim,isign)
        N(1) = Np
        N(2) = Np
        call fourn1(PupilFunc1,N,2,1)
        call fourn1(PupilFunc2,N,2,1)

* The intensity of the Fourier Transform of the Pupil function
* corresponds to an elementary image
*
* Recenter images (necessary for the subsequent convolution ...)
        call jlp_recent_cmplx(PupilFunc1, Np)
        call jlp_recent_cmplx(PupilFunc2, Np)

* Form speckle/fringe patterns
* i.e., intensity of the complex field (PupilFunc1 or PupilFunc2)
        ww = 1. / real(Np * Np)
        do ix=1, Np
          do iy=1, Np
            PSF1(ix,iy) = PSF1(ix,iy) +
     1         ww*real(PupilFunc1(ix,iy)*conjg(PupilFunc1(ix,iy)))
            PSF2(ix,iy) = PSF2(ix,iy) +
     2         ww*real(PupilFunc2(ix,iy)*conjg(PupilFunc2(ix,iy)))
          end do
        end do

  100   continue
******************* END OF LOOP1 (on nsamp) ********************

C Output PSF1:
        if(output_screen.eq.1)then
          out_file='s_psf'
          out_comments='Intensity of complex field in image plane'
          call jlp_writeimag(PSF1,Np,Np,Np,out_file,out_comments)
          output_screen = 2
        endif

        return
        end
********************************************************************
* Compute elementary images by convolving PSF1 with source
* and add the corresponding power spectrum to pwr
*
* INPUT:
* PSF1: elementary PSF
* FT_object: Fourier Transform of the source
* Np: size of power spectrum (=Ns/2)
*
* INPUT/OUTPUT:
* pwr(Np,Np): power spectrum
********************************************************************
       subroutine process_elementary_image(pwr, PSF1,
     1                      FT_object, Np, output_screen)
       implicit none
       integer Np, output_screen
       real*4  pwr(Np,Np), PSF1(Np,Np)
       complex FT_object(Np,Np)

* Temporary arrays
       real*4 ary1(Np,Np)
       complex cary1(Np,Np), TransferFunc(Np,Np)

* Miscellaneous:
       real*4 work(Np,Np)
       integer N(2), jx, jy, iz, jz
       character out_file*60, out_comments*80

        N(1) = Np
        N(2) = Np

C JLP2008
        do jx=1,Np
          do jy=1,Np
             TransferFunc(jx,jy) = complex(PSF1(jx,jy),0.0)
          end do
        end do
        call fourn1(TransferFunc,N,2,1)

C Convolution is done in Fourier domain:
        do jx=1,Np
          do jy=1,Np
C Product is a complex product here!
             cary1(jx,jy) = TransferFunc(jx,jy) * FT_object(jx,jy)
          end do
        end do

        do jx=1,Np
        do jy=1,Np
          pwr(jx,jy) = pwr(jx,jy)
     1                 + real(cary1(jx,jy)*conjg(cary1(jx,jy)))
        end do
        end do

C Output of an elementary frame:
        if(output_screen.eq.2)then

        call fourn1(cary1,N,2,-1)
        do jx=1,Np
          do jy=1,Np
            ary1(jx,jy) = real(cary1(jx,jy))
          end do
        end do
        call jlp_recent_real(ary1, Np)

          do iz=1, Np
            do jz=1, Np
              work(iz, jz) = ary1(iz,jz)
            end do
          end do
          out_file='s_elem_frame'
          out_comments='PSF1 convolved by source'
          call jlp_writeimag(work,Np,Np,Np,out_file,out_comments)
          output_screen = 3
        endif

        return
        end


text:
module parff_parsed_item_m
    use parff_parsed_value_m, only: parsed_value_t

    implicit none
    private
    public :: parsed_item_t

    type :: parsed_item_t
        class(parsed_value_t), allocatable :: item
    end type

    interface parsed_item_t
        module procedure constructor
    end interface
contains
    function constructor(item) result(parsed_item)
        class(parsed_value_t), intent(in) :: item
        type(parsed_item_t) :: parsed_item

        allocate(parsed_item%item, source = item)
    end function
end module

max_stars_repo_name: everythingfunctional/parff

max_issues_repo_name: everythingfunctional/parff

max_forks_repo_name: everythingfunctional/parff


text:
      integer function numbig(x,nx,xmin)
c----------------------------------------------------------------------
c        scans array x for first element less than or equal to
c  cutoff value xmin.  useful primarily for arrays arranged in a
c  decreasing sequence.  returns nx if all elements are larger than
c  xmin as expected.
c
c  arguments-
c    x - array to be scanned
c    nx - length of x
c    xmin - cutoff value to test against.
c----------------------------------------------------------------------
      integer nx
      real x(nx)
      real xmin
      integer i
      do 100 i=1,nx
      if(x(i).lt.xmin) then
            numbig = i - 1
            return
      endif
  100 continue
      numbig = nx
      return
      end

c $Id$

ext: f


text:
/home/kevin/chipyard/sims/verilator/generated-src/chipyard.TestHarness.SmallBoomConfig/EICG_wrapper.v
/home/kevin/chipyard/sims/verilator/generated-src/chipyard.TestHarness.SmallBoomConfig/emulator.cc
-FI /home/kevin/chipyard/sims/verilator/generated-src/chipyard.TestHarness.SmallBoomConfig/verilator.h
/home/kevin/chipyard/sims/verilator/generated-src/chipyard.TestHarness.SmallBoomConfig/SimSerial.cc
/home/kevin/chipyard/sims/verilator/generated-src/chipyard.TestHarness.SmallBoomConfig/testchip_tsi.cc
-FI /home/kevin/chipyard/sims/verilator/generated-src/chipyard.TestHarness.SmallBoomConfig/testchip_tsi.h
/home/kevin/chipyard/sims/verilator/generated-src/chipyard.TestHarness.SmallBoomConfig/SimDRAM.cc
-FI /home/kevin/chipyard/sims/verilator/generated-src/chipyard.TestHarness.SmallBoomConfig/mm.h
/home/kevin/chipyard/sims/verilator/generated-src/chipyard.TestHarness.SmallBoomConfig/mm.cc
-FI /home/kevin/chipyard/sims/verilator/generated-src/chipyard.TestHarness.SmallBoomConfig/mm_dramsim2.h
/home/kevin/chipyard/sims/verilator/generated-src/chipyard.TestHarness.SmallBoomConfig/mm_dramsim2.cc
/home/kevin/chipyard/sims/verilator/generated-src/chipyard.TestHarness.SmallBoomConfig/SimDTM.cc
/home/kevin/chipyard/sims/verilator/generated-src/chipyard.TestHarness.SmallBoomConfig/SimJTAG.cc
-FI /home/kevin/chipyard/sims/verilator/generated-src/chipyard.TestHarness.SmallBoomConfig/remote_bitbang.h
/home/kevin/chipyard/sims/verilator/generated-src/chipyard.TestHarness.SmallBoomConfig/remote_bitbang.cc

max_stars_repo_name: ekiwi/maltese

max_issues_repo_name: ekiwi/maltese

max_forks_repo_name: ekiwi/maltese


text:
module local_routines
    !! Define initial and boundary conditions

    use global_mod, only: dp, ip, charlen, wall_elevation, pi, gravity
    use domain_mod, only: domain_type, STG, UH, VH, ELV
    use read_raster_mod, only: read_gdal_raster
    use which_mod, only: which
    use file_io_mod, only: count_file_lines
    use linear_interpolator_mod, only: linear_interpolator_type
    implicit none

    type :: boundary_information_type
        !! Hold some data used by the boundary condition. We can set this from
        !! inside the main program.
        real(dp) :: offshore_elev
        real(dp) :: boundary_wave_period
    end type

    type(boundary_information_type), public :: boundary_information
        !! The main program will modify this type to set up the boundary condition

    contains

    pure function boundary_function(domain, t, i, j) result(stage_uh_vh_elev)
        !! Make a function to evaluate the boundary at the domain. We will use
        !! this in conjunction with a flather type radiation condition

        type(domain_type), intent(in):: domain
        real(dp), intent(in):: t
        integer(ip), intent(in) :: i, j
        real(dp):: stage_uh_vh_elev(4), wavelength, waveperiod, theta

        ! Specify a wave, corresponding to the far-field condition
        waveperiod = boundary_information%boundary_wave_period
        wavelength = (sqrt(-gravity*boundary_information%offshore_elev)*waveperiod)

        ! Make the wave come from the east
        theta = max((t/waveperiod  + (domain%x(i)-domain%x(domain%nx(1)))/wavelength), 0.0_dp)
        stage_uh_vh_elev(STG) = 1.0_dp * sin(2*pi* theta)
        stage_uh_vh_elev(UH) = -wavelength/waveperiod * stage_uh_vh_elev(STG)
        stage_uh_vh_elev(VH) = 0.0_dp
        stage_uh_vh_elev(ELV) = boundary_information%offshore_elev

    end function

    subroutine set_initial_conditions_circular_island(domain, offshore_depth, island_radius, slope_radius)
        !! Setup initial conditions + locations of gauges

        class(domain_type), target, intent(inout):: domain
        real(dp), intent(in) :: offshore_depth, island_radius, slope_radius

        ! Local variables
        real(dp) :: x, y, elev, slope, radius, theta
        real(dp), allocatable :: gauges(:,:)
        integer(ip) :: i, j

        ! Stage, UH, VH
        domain%U(:,:,[STG, UH, VH]) = 0.0_dp

        ! Set elevation
        slope = offshore_depth/(slope_radius - island_radius)
        do j = 1, domain%nx(2)
            do i = 1, domain%nx(1)
                x = domain%x(i)
                y = domain%y(j)
                radius = sqrt(x*x + y*y)
                elev = -offshore_depth + max(slope_radius - radius, 0.0_dp) *  slope
                domain%U(i,j,ELV) = elev
            end do
        end do
        ! No negative depth!
        domain%U(:,:,STG) = max(domain%U(:,:,STG), domain%U(:,:,ELV))

        ! Set the tide gauges
        allocate(gauges(2, 150))

        ! Gauges around the island (slightly increase the radius to avoid points on dry cells)
        do i = 1, 50
            theta = 2*pi/50.0_dp * (i-1)
            gauges(1:2,i) = (island_radius + domain%dx(1)) *  [cos(theta), sin(theta)]
        end do

        ! Gauges a bit further offshore
        do i = 1, 50
            theta = 2*pi/50.0_dp * (i-1)
            gauges(1:2, i+50) = (island_radius + (slope_radius - island_radius)/2.0_dp) *  [cos(theta), sin(theta)]
        end do

        ! Gauges yet further offshore
        do i = 1, 50
            theta = 2*pi/50.0_dp * (i-1)
            gauges(1:2, i+100) = (slope_radius) *  [cos(theta), sin(theta)]
        end do

        call domain%setup_point_gauges(gauges)

    end subroutine

end module

!@!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

program circular_island
    !! Plane wave scattering around a circular conical island.
    !! This case has an analytical solution, due to:
    !! Zhang and Zhu (1994) New solutions for the propagation of long waves over
    !! variable depth. Journal of fluid mechanics 278: 391-406

    use global_mod, only: ip, dp, minimum_allowed_depth, default_linear_timestepping_method
    use domain_mod, only: domain_type
    use boundary_mod, only: boundary_stage_transmissive_normal_momentum, flather_boundary
    use linear_interpolator_mod, only: linear_interpolator_type
    use local_routines
    implicit none

    integer(ip):: j
    real(dp):: last_write_time, rain_rate
    type(domain_type):: domain

    ! Approx timestep between outputs
    real(dp), parameter :: approximate_writeout_frequency = 50.0_dp
    real(dp), parameter :: final_time = 200000.0_dp

    ! Domain info
    character(charlen) :: timestepping_method = default_linear_timestepping_method

    !! length/width
    real(dp), dimension(2) :: global_lw, global_ll
    integer(ip), dimension(2) :: global_nx

    ! Local variables
    real(dp) :: timestep
    real(dp) :: offshore_depth, island_radius, slope_radius, dx
    real(dp) :: boundary_wave_period
    character(charlen):: test_case, bc_file
    real(dp):: tank_bases(4), tank_slopes(4)
    integer(ip) :: full_write_step

    ! Write the stage raster time-series less often than we write at gauges, to avoid
    ! overly large files.
    integer(ip), parameter :: frequency_full_write_steps = 60
    logical, parameter :: never_write_grid_time_slices = .false.

    ! Zero the max stage record after this much time has elapsed. The
    ! idea is to allow the transients to pass, then reset the max stage,
    ! so it ultimately only records the 'stationary' part of the run
    real(dp), parameter :: reset_max_stage_at_time = 120000.0_dp

    ! Resolution
    dx = 2000.00_dp

    ! Large scale domain. Wave comes from east side
    global_lw = [2000.0_dp , 3000.0_dp]*1e+03
    global_ll = -global_lw/2.0_dp
    global_nx = nint(global_lw/dx)

    ! Geometry
    island_radius = 40000.0_dp
    slope_radius = 160000.0_dp
    offshore_depth = 300.0_dp

    ! Incoming wave
    ! wavelength = 2 x slope_radius
    boundary_wave_period = slope_radius * 2.0_dp /sqrt(9.8_dp * offshore_depth) !12.0_dp * 60.0_dp

    domain%timestepping_method = timestepping_method

    ! Allocate domain -- must have set timestepping method BEFORE this
    call domain%allocate_quantities(global_lw, global_nx, global_ll)

    ! Call local routine to set initial conditions
    call set_initial_conditions_circular_island(domain, offshore_depth, island_radius, slope_radius)

    ! Get the boundary data and make an interpolation function f(t) for gauge 4
    domain%boundary_function => boundary_function
    domain%boundary_subroutine => flather_boundary
    boundary_information%offshore_elev = -offshore_depth
    boundary_information%boundary_wave_period = boundary_wave_period

    ! Linear requires a fixed timestep
    if (.not. domain%adaptive_timestepping) then
        timestep = domain%stationary_timestep_max()
    end if

    ! Trick to get the code to write out just after the first timestep
    last_write_time = domain%time - approximate_writeout_frequency

    full_write_step = 0
    ! Evolve the code
    do while (.true.)

        if(domain%time - last_write_time >= approximate_writeout_frequency) then

            ! Reset the peak stage record once during the simulation, so the final
            ! result reflects the 'stationary' model solution
            if((domain%time >= reset_max_stage_at_time) .and. &
               (last_write_time <= reset_max_stage_at_time)) then
                domain%max_U(:,:,1) = -huge(1.0_dp)
            end if

            last_write_time = last_write_time + approximate_writeout_frequency

            ! This avoids any artefacts in the numerical update of the model
            ! which should be overwritten by the boundary condition
            call domain%update_boundary()

            full_write_step = full_write_step + 1
            if(mod(full_write_step, frequency_full_write_steps) == 0) then
                call domain%write_to_output_files(time_only=never_write_grid_time_slices)
                call domain%print()
                print*, 'Mass balance: ', domain%mass_balance_interior()
            end if

            call domain%write_gauge_time_series()

        end if

        if (domain%time > final_time) exit

        if (.not. domain%adaptive_timestepping) then
            call domain%evolve_one_step(timestep = timestep)
        else
            call domain%evolve_one_step()
        end if

    end do

    call domain%write_max_quantities()

    ! Print timing info
    call domain%timer%print()

    call domain%finalise()

end program

max_stars_repo_name: KayaGeoAus/ptha

max_issues_repo_name: KayaGeoAus/ptha

max_forks_repo_name: KayaGeoAus/ptha


text:
      subroutine besjn (x, jl, nl)

c-----------------------------------------------------------------------
c
c     purpose:  to calculate the spherical bessel functions jl and nl
c               for l = 0 to 30 (no offset)
c
c     arguments:
c       x = argument of jl and nl
c       jl = jl bessel function (abramowitz conventions)
c       nl = nl bessel function (abramowitz yl conventions)
c            Note that this array nl = abramowitz yl.
c       jl and nl must be dimensioned
c            complex*16 jl(ltot+2), nl(ltot+2), with ltot defined in
c            dim.h.
c
c     notes:  jl and nl should be calculated at least to 10 place
c             accuracy for the range 0<x<100 according to spot
c             checks with tables
c
c     error messages written with PRINT statement.
c
c     first coded by r. c. albers on 14 dec 82
c
c     version 3
c
c     last modified: 27 jan 83 by r. c. albers
c     dimension of jl,nl changed from 31 to 26  (10 aug 89) j. rehr
c     modified again, siz, June 1992
c
c-----------------------------------------------------------------------

       implicit double precision (a-h, o-z)
       include 'dim.h'

       complex*16 x
       complex*16 jl(ltot+2), nl(ltot+2)
       complex*16 cjl(ltot+2), sjl(ltot+2), cnl(ltot+2), snl(ltot+2)

       complex*16 xjl,xnl,asx,acx
       complex*16 xi,xi2,xi3,xi4,xi5,xi6,xi7,xi8,xi9,xi10,xi11
       double precision xcut, xcut1, xcut2

       parameter (xcut = 1.d0, xcut1 = 7.51d0, xcut2 = 5.01d0)

      if (dble(x) .le. 0)  call fstop(' at besjn: '//
     $      'Re(x) is <=0')

      lmaxp1 = ltot+2

      if (dble(x) .lt. xcut)  then
c        case Re(x) < 1, just use series expansion
         do 10 il = 1,lmaxp1
            l = il-1
            ifl = 0
            call bjnser (x,l,xjl,xnl,ifl)
            jl(il) = xjl
            nl(il) = xnl
   10    continue

      elseif (dble(x) .lt. xcut1)  then

c        case 1 <= Re(x) < 7.5

         call bjnser (x,lmaxp1-1,xjl,xnl,1)
         jl(lmaxp1) = xjl

         call bjnser (x,lmaxp1-2,xjl,xnl,1)
         jl(lmaxp1-1) = xjl

         if (dble(x) .lt. xcut2)  then
c           Re(x) < 5
            call bjnser (x,0,xjl,xnl,2)
            nl(1) = xnl
            call bjnser (x,1,xjl,xnl,2)
            nl(2) = xnl
         else
c           Re(x) >= 5
            asx = sin(x)
            acx = cos(x)
            xi = 1 / x
            xi2 = xi*xi
            nl(1) = -acx*xi
            nl(2) = -acx*xi2 - asx*xi
         endif

c        Use recursion relation 10.1.19 to get nl and jl
         do 50 lp1 = 3, lmaxp1
            tlxp1 = 2*lp1 -3
            nl(lp1) = tlxp1 * nl(lp1-1) / x  -  nl(lp1-2)
   50    continue

         do 60 lx = 3,lmaxp1
            lp1 = lmaxp1+1-lx
            l = lp1-1
            tlxp3 = 2*l + 3
            jl(lp1) = tlxp3 * jl(lp1+1) / x  -  jl(lp1+2)
   60    continue

      else
c        case Re(x) > 7.5
c        Use AS 10.1.8 and 10.1.9, sjl=P, qjl=Q, note that AS formulae
c        use cos (z - n*pi/2), etc., so cos and sin terms get a bit
c        scrambled (mod 4) here, since n is integer.  These are hard-
c        coded into the terms below.
         xi = 1 / x
         xi2  = xi*xi
         xi3  = xi*xi2
         xi4  = xi*xi3
         xi5  = xi*xi4
         xi6  = xi*xi5
         xi7  = xi*xi6
         xi8  = xi*xi7
         xi9  = xi*xi8
         xi10 = xi*xi9
         xi11 = xi*xi10

         sjl(1) = xi
         sjl(2) = xi2
         sjl(3) = 3.d0*xi3 - xi
         sjl(4) = 15.d0*xi4 - 6.d0*xi2
         sjl(5) = 105.d0*xi5 - 45.d0*xi3 + xi
         sjl(6) = 945.d0*xi6 - 420.d0*xi4 + 15.d0*xi2
         sjl(7) = 10395.d0*xi7 - 4725.d0*xi5 + 210.d0*xi3 - xi
         sjl(8) = 135135.d0*xi8 - 62370.d0*xi6 + 3150.d0*xi4 -
     $        28.d0*xi2
         sjl(9) = 2027025.d0*xi9 - 945945.d0*xi7 + 51975.d0*xi5 -
     1        630.d0*xi3 + xi
         sjl(10) = 34459425.d0*xi10 - 16216200.d0*xi8 +
     1        945945.d0*xi6 - 13860.d0*xi4 + 45.d0*xi2
         sjl(11) = 654729075.d0*xi11 - 310134825.d0*xi9 +
     1        18918900.d0*xi7 - 315315.d0*xi5 + 1485.d0*xi3 - xi

         cjl(1) = 0
         cjl(2) = -xi
         cjl(3) = -3.d0*xi2
         cjl(4) = -15.d0*xi3 + xi
         cjl(5) = -105.d0*xi4 + 10.d0*xi2
         cjl(6) = -945.d0*xi5 + 105.d0*xi3 - xi
         cjl(7) = -10395.d0*xi6 + 1260.d0*xi4 - 21.d0*xi2
         cjl(8) = -135135.d0*xi7 + 17325.d0*xi5 - 378.d0*xi3 + xi
         cjl(9) = -2027025.d0*xi8 + 270270.d0*xi6 - 6930.d0*xi4 +
     $        36.d0*xi2
         cjl(10) = -34459425.d0*xi9 + 4729725.d0*xi7 - 135135.d0*xi5 +
     1        990.d0*xi3 - xi
         cjl(11) = -654729075.d0*xi10 + 91891800.d0*xi8 -
     1        2837835.d0*xi6 + 25740.d0*xi4 - 55.d0*xi2
         do 80 ie = 1,11
            snl(ie) = cjl(ie)
            cnl(ie) = -sjl(ie)
   80    continue
         do 90 lp1 = 12,lmaxp1
            tlxp1    = 2*lp1-3
            sjl(lp1) = tlxp1*xi*sjl(lp1-1)-sjl(lp1-2)
            cjl(lp1) = tlxp1*xi*cjl(lp1-1)-cjl(lp1-2)
            snl(lp1) = tlxp1*xi*snl(lp1-1)-snl(lp1-2)
            cnl(lp1) = tlxp1*xi*cnl(lp1-1)-cnl(lp1-2)
   90    continue
         asx = sin(x)
         acx = cos(x)
         do 110 lp1 = 1,lmaxp1
            jl(lp1) = asx*sjl(lp1)+acx*cjl(lp1)
            nl(lp1) = asx*snl(lp1)+acx*cnl(lp1)
  110    continue
      endif

      return
      end

max_stars_repo_name: keechul/ifeffit

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_name: bruceravel/ifeffit

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_name: bruceravel/ifeffit

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


text:
C
C $Id: slgeti.f,v 1.6 2008-07-27 00:17:26 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE SLGETI (PNAM,IVAL)
C
C Get in IVAL the integer value of the STITLE parameter named PNAM.
C
        CHARACTER*(*) PNAM
C
C Check for an uncleared prior error.
C
        IF (ICFELL('SLGETI - UNCLEARED PRIOR ERROR',1).NE.0) RETURN
C
C Call SLGETR to obtain the real value of the parameter and then
C return the integer portion of that.
C
        CALL SLGETR (PNAM,RVAL)
        IF (ICFELL('SLGETI',2).NE.0) RETURN
        IVAL=INT(RVAL)
C
C Done.
C
        RETURN
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


!        Generated by TAPENADE     (INRIA, Ecuador team)
!  tapenade 3.x
!
!  Differentiation of omp_old_save in reverse (adjoint) mode (with options noISIZE i8):
!   gradient     of useful results: *concold *temp *headold *head
!                *tempold *presold *conc *pres
!   with respect to varying inputs: *concold *temp *headold *head
!                *tempold *presold *conc *pres
!   Plus diff mem management of: concold:in temp:in headold:in
!                head:in tempold:in presold:in conc:in pres:in
!> @brief save current state as an old version
!> @param[in] level level number (which old version)
!> @param[in] ismpl local sample index
SUBROUTINE OMP_OLD_SAVE_AD(level, ismpl)
  use arrays

  USE ARRAYS_AD

  USE MOD_GENRL
  USE MOD_CONC
  IMPLICIT NONE
! local sample index
  INTEGER :: ismpl
! cgen level index
  INTEGER, INTENT(IN) :: level
! directional cell-indices
  INTEGER :: i, j, k
! counter for concentration
  INTEGER :: l
! Start position in array for process
  INTEGER :: tpos
! Number of array elements for process
  INTEGER :: tanz
  INTRINSIC ABS
  INTRINSIC MAX
  INTEGER :: abs0
  INTEGER :: abs1
  INTEGER :: abs2
  INTEGER :: abs3
  INTEGER :: max1
  INTEGER :: max2
  INTEGER :: max3
  INTEGER :: max4
  INTEGER :: branch
! OpenMP partition, get tpos, tanz
  CALL OMP_PART(i0*j0*k0, tpos, tanz)
! Get i, j, k indices for tpos
  CALL IJK_M(tpos, i, j, k)
  IF (ismpl .GE. 0.) THEN
    abs0 = ismpl
  ELSE
    abs0 = -ismpl
  END IF
  IF (1 .LT. ismpl) THEN
    max1 = ismpl
  ELSE
    max1 = 1
  END IF
! Copy all variable arrays to old: var -> varold
  IF (ismpl .GE. 0.) THEN
    abs1 = ismpl
  ELSE
    abs1 = -ismpl
  END IF
  IF (1 .LT. ismpl) THEN
    max2 = ismpl
  ELSE
    max2 = 1
  END IF
  IF (ismpl .GE. 0.) THEN
    abs2 = ismpl
  ELSE
    abs2 = -ismpl
  END IF
  IF (1 .LT. ismpl) THEN
    max3 = ismpl
  ELSE
    max3 = 1
  END IF
  DO l=1,ntrans
    IF (ismpl .GE. 0.) THEN
      CALL PUSHINTEGER8(abs3)
      abs3 = ismpl
      CALL PUSHCONTROL1B(1)
    ELSE
      CALL PUSHINTEGER8(abs3)
      abs3 = -ismpl
      CALL PUSHCONTROL1B(0)
    END IF
    IF (1 .LT. ismpl) THEN
      CALL PUSHINTEGER8(max4)
      max4 = ismpl
      CALL PUSHCONTROL1B(0)
    ELSE
      CALL PUSHINTEGER8(max4)
      max4 = 1
      CALL PUSHCONTROL1B(1)
    END IF
  END DO
  DO l=ntrans,1,-1
    CALL DCOPY_AD(tanz, conc(i, j, k, l, abs3), conc_ad(i, j, k, l, abs3&
&           ), 1, concold(tpos, l, level, max4), concold_ad(tpos, l, &
&           level, max4), 1)
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 0) THEN
      CALL POPINTEGER8(max4)
    ELSE
      CALL POPINTEGER8(max4)
    END IF
    CALL POPCONTROL1B(branch)
    IF (branch .EQ. 0) THEN
      CALL POPINTEGER8(abs3)
    ELSE
      CALL POPINTEGER8(abs3)
    END IF
  END DO
  CALL DCOPY_AD(tanz, pres(i, j, k, abs2), pres_ad(i, j, k, abs2), 1, &
&         presold(tpos, level, max3), presold_ad(tpos, level, max3), 1)
  CALL DCOPY_AD(tanz, temp(i, j, k, abs1), temp_ad(i, j, k, abs1), 1, &
&         tempold(tpos, level, max2), tempold_ad(tpos, level, max2), 1)
  CALL DCOPY_AD(tanz, head(i, j, k, abs0), head_ad(i, j, k, abs0), 1, &
&         headold(tpos, level, max1), headold_ad(tpos, level, max1), 1)
END SUBROUTINE OMP_OLD_SAVE_AD


text:
\ SPI flash access

\ Flash addresses are doubles to support 16-bit and 18-bit cell size.
\ Compiled strings (for `type` etc.) are assumed to have 0 for the upper cell
\ of the address. Those strings must be below flash address 2^cellsize.

there
hex

: _isp  ( c -- )                \ write ISP byte to SPIF
   io'isp io!
;
: ispwait  ( -- )               \ wait for ISP command to finish
   begin  io'isp io@  while  noop  repeat
;
: ispcmd  ( c -- )              \ c --
   _isp ispwait                 \ write ISP command
;
: ispnum  ( n -- )
   3F and _isp                  \ 6-bit number command
;

\ The gecko key is loaded by writing to a register and then triggering a load.
\ `fcmd24` is used for writing to flash, which isn't done currently

: gkey                          \ n --
   io'gkey io!                  \ shift in a new key
;
: )gkey  ( -- )                 \ --
   48 ispcmd                    \ re-key the gecko keystream
;
: fcmd24  ( df-addr cmd len --) \ set start address and set command
   0 _isp _isp  82 _isp  ispcmd \ len is 3 + extra chars ( df-addr )
   2dup 10 drshift drop
   BASEBLOCK +   ispcmd         \ addr[23:16]
   over swapb ispcmd            \ addr[15:8]
   drop ispcmd                  \ addr[7:0]
;

: )@f  ( -- )                   \ 2.4030 --
   80 _isp                      \ end read (raise CS#)
;

\ Load a boot stream from flash starting at a specified page.
\ If the data is glitched, such as from an ESD hit, the CRC32 won't match
\ so retry until it's okay.

: spifload  ( page -- )
   api !
   begin
      api @
      dup 6 rshift ispnum ispnum \ set page number
      44 ispcmd                  \ interpret the flash stream
   io'boot io@ bootokay and until
;

\ Large applications are supported by caching, which is enabled by an exception.
\ The `api_recover` exception is triggered when R>1FFFh.
\ Instead of popping PC from the return stack, execution jumps to a fixed addr.
\ `api_recover` gets the required page from R and masks off the page number.

\ `xexec` executes a word at a specified xxt (extended execution token)
\ consisting of the flash page number packed with a 13-bit address.
\ in the format 11.13 where the upper 11 bits are the flash page number
\ (a page is 256 bytes) and the lower 13 bits are the xt.

: xexec  ( xxt -- ) ( R: addr -- page|addr )
   api @ 0d lshift r> + >r              \ add page to return address
   1fff invert overand                  ( xxt new_page )
   0d rshift spifload
   1fff and  >r                         \ execute
; no-tail-recursion
' xexec resolves api_trap

\ still has an occasional problem with return stack underflow
\ `xcall` saves the return address into old code, which `api_trap` appends a
\ current page number to. The new page is loaded and execution starts at the
\ new page's code. Upon execute, the return stack has the old page number.

\ When a return address with page<>0 is encountered, the pc jumps to api_recover
\ while the questionable return address is on the return stack.
\ The old page is reloaded before execution returns.

:noname  ( R: page|addr -- addr )       \ returning to an api word (page <> 0)
   r> 1fff overand >r
   0d rshift  spifload                          ( needed_page )
\   dup api @ xor if                     \ returning to a corrupted page
\      spifload exit
\   then drop
; resolves api_recover

decimal

\ Flash reads using hardware instead of discrete SPI transfers.
\ Double cell addresses are used but not implemented. If they are needed,
\ you can add code to handle the upper half of the address.

: fwait  ( -- )                 \ wait for flash read to finish
   begin  io'gkey io@  while  noop  repeat
;
: _x@f  ( df-addr cfg -- x )    \ hardware flash read
   nip  io'fcfg io!             \ flash read setup
   io'fread io!                 \ trigger a read
   fwait  io'fread io@          \ get result
;
: c@f>  ( -- c )
   [ 2 2* 2* ] literal
   io'fcfg io!                  \ flash read setup
[ ;
: @f>  ( -- x )                 \ read using the default setup
   io'fnext dup io!             \ trigger a read-next
   fwait  io'fread io@          \ get result
;
: @f(  ( df-addr -- x )         \ 3-byte big-endian read
   [ 2 2* 2*  2 + ] literal     \ format=2, size=2
   _x@f
;
: @f  ( df-addr -- x )          \ 3-byte big-endian read
   @f(  )@f                     \ raise CS afterwards
;
: w@f(  ( df-addr -- x )        \ first 2-byte big-endian read
   [ 2 2* 2*  1 + ] literal
   _x@f
;
: w@f  ( df-addr -- x )         \ 2-byte big-endian read
   w@f(  )@f
;
: c@f(  ( df-addr -- x )        \ first 1-byte read
   [ 2 2* 2* ] literal
   _x@f
;
: c@f  ( df-addr -- x )         \ 1-byte read
   c@f(  )@f
;

: d@f    @f 0 ;                 \ 2.4055 df-addr -- d

: fcount                        \ 2.4060 df-addr -- df-addr+1 c
   2dup 1 0 d+  2swap c@f
;
: 3*   ( n -- 3n ) dup 2* + ;   \ multiply by 3, goes with f@

\ `emit` may require reading font bitmaps from flash, so `f$type` reads the
\ string into a RAM buffer because the bitmap read disrupts keystream sync.
\ This is a good place for `pad`. Make it big enough for app usage.
\
\ `fbuf` moves a string from flash to RAM.
\ The keystream is assumed to be in sync.

: fbuf                         \ 2.4065 df-addr c-addr u --
   >r >r  c@f(  r@ c!  r> r>   ( c-addr u ) \ 1st char
   begin  1 /string  dup while
      over @f> swap c!
   repeat  2drop  )@f
;

256 equ |pad|
|pad| buffer: pad

: ftype                         \ 2.4070 df-addr u --
   >r pad r@ fbuf  pad r> type
;

tkey or [if]
: /text  ( f-addr -- df-addr )  \ synchronize keystream
   tkey gkey gkey  dup gkey )gkey  0
;
[else]
: /text  ( f-addr -- df-addr )
   0
;
[then]

: f$type                        \ 2.4080 f-addr --
   /text  fcount  ftype         \ emit the "flash string"
;

: [wid]   ( wid -- addr )               \ wid is indexed from 1
   cells [ wids cell - ] literal +
;
: wid    [wid] @ ;                      \ wid -- f-addr

\ | Length  | Name  | Usage                |
\ | ------- |:-----:| --------------------:|
\ | M bytes | link  | Link to next header  | <-- keystream
\ | 1 byte  | N     | Length of `name`     |
\ | N bytes | name  | Name string          |
\ | M bytes | xte   | *xt* for Execution   | <-- ht
\ | M bytes | xtc   | *xt* for Compilation |
\ | M bytes | w     | Optional data        |
\ | 1 byte  | flags | Packed flags (0=on)  |
\ | 1 byte  | app   | Applet ID            |
\
\ For `find`, the SPI read sequence is continued until it doesn't match.
\ `match` compares a key in data space to the flash byte stream.
\ It's a case-sensitive comparison. For a case-insensitive search, insert a
\ tolower function.

: match  \ addr len -- found?
   for                                  \ addr' | len'
      count ( tolower ) @f> xor if
         rdrop dup xor exit             \ mismatch in string
      then
   next  drop 1                         \ all bytes have been read
;

\ Search one wordlist, returning the address of header data.
: _hfind  \ 2.4100 addr len wid -- addr len 0 | addr len ht
   wid  over 31 u> -19 and throw        \ name too long
   over 0= if dup xor exit then         \ addr 0 0   zero length string
   begin
      dup>r  /text
      @f( >r  c@f>                      \ addr len1 len2 | head link
      over = if
         2dup match if                  \ found
            rdrop dup r> +              \ 'link + NameLength + 1 + cellbytes
            [ cellbits 7 + 8 /  1 + ] literal +
            [ 2 2* 2*  2 + ] literal    \ next reads will be 24-bit
            io'fcfg io!   exit          \ don't end the read yet
         then
      then )@f r>  rdrop                \ end flash read
   dup 0= until                         \ not found
;

\ A primitive for `find` that returns a ht.
: hfind  \ 2.4110 addr len -- addr len | 0 ht  \ search the search order
   #order @  begin  dup  while  1-      \ addr len idx
      dup>r  cells orders + @ _hfind    \ addr len 0/ht
      ?dup if
         rdrop 2drop 0 swap exit        \ found, exit
      then
   r> repeat  drop
;

there swap - . .( instructions used by flash) cr

max_stars_repo_name: bradleyeckert/chad

max_issues_repo_name: bradleyeckert/chad

max_forks_repo_name: bradleyeckert/chad


text: While it used to be one of the modes of transportation in town, it is now a Departed Businesses departed business. They offered special rates to the airport.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search