!some input/output format for the theta code:
     ! set initial conditions of the synchronous particle:
      open(unit=13,file='sim.dat',status='old')
      read(13,*)frqhz,ekk,pmass,a1in,charge,idata,ierr,Iphase,nsizeperd, &
         refphin,tolin,deltaxin,refphinB,tolinB,deltaxinB
..... 
        read(13,*,end=22)zlen,ns,mapstp,itype,val1,val2,val3,val4,&
                          val5,val6,val7,val8,val9,val10,val11
.....
        print*,"Nquad, Nrf: ",Nquad,Nrf
        print*,"input random seed:"
        read(*,*)nrandseed
        !seedarray(1) = 1001+5*nrandseed
        !call random_seed(PUT=seedarray(1:1))
        call random_number(testrand)
        print*,"testrand",testrand
        open(unit=14,file='errorlimit.dat',status='old')
        !read in the limits for quadrupole errors
        read(14,*)qxlmt,qylmt,qzlmt,qrotxlmt,qrotylmt,qrotzlmt,qfldlmt
        tmpsum1 = 0.0
        tmpsum2 = 0.0
        tmpsum3 = 0.0
        tmpsum4 = 0.0
        tmpsum5 = 0.0
        tmpsum6 = 0.0
        tmpsum7 = 0.0
        tmpsum8 = 0.0
        do i = 1, Nrf
          call random_number(rd)
          rfxerr(i) = rfxlmt*(2*rd - 1.0)
          call random_number(rd)
          rfyerr(i) = rfylmt*(2*rd - 1.0)
          call random_number(rd)
          rfzerr(i) = rfzlmt*(2*rd - 1.0)
          call random_number(rd)
          rfrotxerr(i) = rfrotxlmt*(2*rd - 1.0)
          call random_number(rd)
          rfrotyerr(i) = rfrotylmt*(2*rd - 1.0)
          call random_number(rd)
          rfrotzerr(i) = rfrotzlmt*(2*rd - 1.0)
          call random_number(rd)
          rfflderr(i) = rffldlmt*(2*rd-1.0)
          call random_number(rd)
          print*,"rd: ",i,rd,rfphlmt
          rfpherr(i) = rfphlmt*(2*rd-1.0)
          tmpsum1 = tmpsum1 + rfxerr(i)
          tmpsum2 = tmpsum2 + rfyerr(i)
          tmpsum3 = tmpsum3 + rfzerr(i)
          tmpsum4 = tmpsum4 + rfflderr(i)
          tmpsum5 = tmpsum5 + rfpherr(i)
          tmpsum6 = tmpsum6 + rfrotxerr(i)
          tmpsum7 = tmpsum7 + rfrotyerr(i)
          tmpsum8 = tmpsum8 + rfrotzerr(i)
        enddo
        do i = 1, Nrf
          rfxerr(i) = rfxerr(i) - tmpsum1/Nrf
          rfyerr(i) = rfyerr(i) - tmpsum2/Nrf
          rfzerr(i) = rfzerr(i) - tmpsum3/Nrf
          rfflderr(i) = rfflderr(i) - tmpsum4/Nrf
          rfpherr(i) = rfpherr(i) - tmpsum5/Nrf
          rfrotxerr(i) = rfrotxerr(i) - tmpsum6/Nrf
          rfrotyerr(i) = rfrotyerr(i) - tmpsum7/Nrf
          rfrotzerr(i) = rfrotzerr(i) - tmpsum8/Nrf
        enddo
        close(14)
        rewind(13)
        read(13,*)frqhz,ekk,pmass,a1in,charge,idata,ierr,Iphase,nsizeperd, &
                  refphin,tolin,deltaxin,refphinB,tolinB,deltaxinB
      endif
 
      open(3,file="fort.18",status="unknown")
101   continue
      read(3,*,end=111)val1,val2,val3,val4,val5,val6
      a1in = val2
      ekk = val4*1.0e6
      goto 101
111   continue
      close(3)
      print*,"a1in, ekk: ",a1in,ekk
.....
      !initialize beam line element information:
      !----------------------------------------------
      read(13,*,end=999)zlen,ns,mapstp,itype,val1,val2,val3,val4,&
                        val5,val6,val7,val8,val9,val10,val11
......
!-------------------------------------------------------------
!<<< The following is used for reading discrete data
!-------------------------------------------------------------
        if(idata.eq.1) then  !use discrete data
          read(8,*,end=999)arg1,arg2,arg3
          n=n+1
          zdat(n)=arg1
          edat(n)=arg2
          epdat(n)=arg3
          if(Emax.le.abs(arg2)) then
            Emax = abs(arg2)
          endif
!-------------------------------------------------------------
!<<< The following is used for reading Fourier coefficients.
!-------------------------------------------------------------
        else !use analytical function
          read(8,*,end=999)arg1
          n=n+1