text:
C                 ******************************
C                      I L L U S T R A T E
C                   Biomolecular Illustration
C		  ******************************
C		  copyright 2019 David S Goodsell
C
C Licensed under the Apache License, Version 2.0 (the "License");
C you may not use this file except in compliance with the License.
C You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0
C Unless required by applicable law or agreed to in writing, software
C distributed under the License is distributed on an "AS IS" BASIS,
C WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
C See the License for the specific language governing permissions and
C limitations under the License
C
C This work was supported by Damon Runyon-Walter Winchell Cancer Research Fund Fellowship DRG 972,
C the US National Institutes of Health R01-GM120604 and the kind support of the RCSB Protein Data Bank.
C
C--------------------------------------------------------------------
C
C DS Goodsell & AJ Olson (1992) "Molecular Illustration in Black and White" JMolGraphics 10, 235-240.
C April 2019 -- Simplified and released with only non-photorealistic rendering
C
C--------------------------------------------------------------------
C compile: gfortran illustrate.f -o illustrate
C run: illustrate < command_file
C
C Command file is read from unit 5
C Command file has command cards, followed by parameter cards
C Idiosyncracies (warning, postdoc code):
C *** must issue command cards in this order
C *** any number of rotation cards may be added, and they are concatenated when added
C      this means they are effectively applied last to first
C      so if you're progressively refining a position, add new rotations to the top of the list
C *** origin at upper left, +x down, +y left to right, +z towards viewer, molecules clipped at z=0
C
C read                                          #READ command
C 2hhb.pdb                                      #PDB format coordinate file
C HETATM-----HOH-- 0,9999, 0.5,0.5,0.5, 0.0     #selection/rendering cards
C ATOM  -H-------- 0,9999, 0.5,0.5,0.5, 0.0
C ATOM  H--------- 0,9999, 0.5,0.5,0.5, 0.0
C ATOM  -C-------A 0,9999, 1.0,0.6,0.6, 1.6
C ATOM  -S-------A 0,9999, 1.0,0.5,0.5, 1.8
C ATOM  ---------A 0,9999, 1.0,0.5,0.5, 1.5
C ATOM  -C-------C 0,9999, 1.0,0.6,0.6, 1.6
C ATOM  -S-------C 0,9999, 1.0,0.5,0.5, 1.8
C ATOM  ---------C 0,9999, 1.0,0.5,0.5, 1.5
C ATOM  -C-------- 0,9999, 1.0,0.8,0.6, 1.6
C ATOM  -S-------- 0,9999, 1.0,0.7,0.5, 1.8
C ATOM  ---------- 0,9999, 1.0,0.7,0.5, 1.5
C HETATMFE---HEM-- 0,9999, 1.0,0.8,0.0, 1.8
C HETATM-C---HEM-- 0,9999, 1.0,0.3,0.3, 1.6
C HETATM-----HEM-- 0,9999, 1.0,0.1,0.1, 1.5
C END                                          #end of READ commad
C center                                       #CENTER command
C auto                                         #use autocentering (typical)
C trans                                        #TRANSLATION command
C 0.,0.,0.                                     #x,y,z for translation
C scale                                        #SCALE command
C 12.0                                         #scale value (pixels/Angstrom)
C zrot                                         #ROTATION command
C 90.                                          #rotation angle (deg)
C wor                                          #WORLD rendering parameter
C 1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0              #rgb background, rgb fog, fractional opacity of fog front & back
C 1,0.0023,2.0,1.0,0.2                         #soft shadow parameters
C -30,-30                                      #image size, negative for autosizing with padding
C illustrate                                   #ILLUSTRATION command
C 3.0,10.0,4,0.0,5.0                           #contour outlines
C 3.0,10.0                                     #subunit outlines
C 3.0,8.0,6000.                                #residue outlines
C calculate                                    #CALCULATE command
C 2hhb.pnm                                     #image file name in ppm format
C
C--------------------------------------------------------------------
c ***** OUTPUT SCAN LINE *****
	integer*4 scanline(9000)
c ***** FRAME BUFFER WITH COLORS *****
	real*4 pix(-10:3008,-10:3008,4)
c ***** FRAME BUFFER CONTAINING THE z HEIGHT *****
	real*4 zpix(-10:3008,-10:3008)
c ***** FRAME BUFFER CONTAINING THE ATOM NUMBER *****
	integer*4 atom(-10:3008,-10:3008)
	integer*4 bio(-10:3008,-10:3008)
	integer*4 n,ia
c ***** SHADING MAPS FOR THE EIGHT ATOM TYPES *****
	real*4 sphdat(0:32000,3)
	integer*4 numpix
c ***** PHONG SHADING PARAMETERS ******
	real*4 colortype(0:1000,3)
	real*4 rback(3),rfog(3)
c ***** ATOMIC INFORMATION *****
	real*4 coord(3,350000)
	integer*4 type(350000),res(0:350000),su(0:350000)
	real*4 radtype(1000)
	character chain,chainlast
c ***** transformation matrices *****
	real*4 matrixin(4,4),rm(4,4)
c ***** biological unit stuff *****
	real*4 biomat(4,4,500)
	integer*4 nbiochain,nbiomat
	character biochain(500)
c ***** STUFF FOR OUTLINES *****
	real*4 l_opacity,l_opacity_ave,g_opacity,opacity
	real*4 l(-1:1,-1:1),g
	real*4 l_low,l_high,g_low,g_high
	real*4 l_diff_min, l_diff_max
c ***** Conical shadows *****
	real*4 rtable(-51:51,-51:51),coneangle,pcone,rcone
c ***** ETC. *****
	real*4 x,y,z,rx,ry,rz,xp,yp,zp,d
	real*4 xn,yn,zn,ci,xs,xy,xz,cs
	character*3 comcode(10),command
	character*80 filename,inputfile
	integer*4 ixsize,iysize, idepth
	integer*4 ix,iy,iz
	integer illustrationflag
c ***** stuff for reading atoms
	integer*4 resrange(2,1000),inptype(1000),inpsu(1000)
	character*6 atomdescriptor(1000)
	character*10 descriptor(1000)
	character*80 instring
c--------------------------------------------------------------------
c ***** Initialize a few things *****
c--------------------------------------------------------------------
	rscale=1.

	illustrationflag=0

	data comcode/'rea','tra','xro','yro','zro','sca','cen',
     &		             'wor','cal','ill'/

	call clearmatrix(rm)
	l_low=1.
	l_high=10.
	g_low=350000.
	g_high=21000.
	xshmax=4000
	yshmax=4000
	su(0)=9999
	res(0)=9999
	xtran=0.
	ytran=0.
	ztran=0.
	l_diff_max=50.
	l_diff_min=1.
	ixsize=0
	iysize=0
	do i=1,3
	do j=1,4
	do k=1,500
	biomat(i,j,k)=0.
	if (i.eq.j) biomat(i,j,k)=1.
	enddo
	enddo
	enddo
c ********************* READ CONTROL CARDS ***************************
 10	read (5,101,end=999) command
	icommand=10
	do icount=1,10
 100	if (command.eq.comcode(icount)) icommand=icount
	enddo
  	goto (1,2,3,4,5,6,7,8,111,12),icommand
	write (6,102) ' ***** invalid control card read: ',
     &       command, ' ***** '
	goto 10
 101	format (a3)
 102	format (a35,a3,a7)
c--------------------------------------------------------------------
c  *** read and classify atoms ***
c
 1	continue
	read(5,113) inputfile
	open(1,file=inputfile,form='formatted',status='old')
c --- read atom descriptors ---
c param: atom descriptor cards
c ATOM  -CA--ALA-A 0,9999,0.5,0.5,0.5,1.6
c ^^^^^^   ATOM or HETATM, matched with columns 1-6 of PDB record
c       ^^^^^^^^^^ compared to columns 13-22 of PDB record
c                  "-" is a wildcard
c                  ^^^^^^ residue number range
c                         ^^^^^^^^^^^ r,g,b (0.-1.)
c                                     ^^^ radius (Angstrom)
c cards are read in order, and if there is a match, the atom is assigned that type
c cards are read until one is found without "ATOM  " or "HETATM"
c
	ndes=0
c --- default 50% gray ---
	do i=0,1000
	do j=1,3
	colortype(i,j)=.5
	enddo
	enddo
 7020	read(5,7100) instring
 	if (instring(1:3).eq.'END') goto 7030
 	  ndes=ndes+1
 	  read(instring,21) atomdescriptor(ndes),descriptor(ndes)
 	  read(instring(18:80),*) (resrange(i,ndes),i=1,2),
     &                     rr,rg,rb,rad
	  colortype(ndes,1)=rr
	  colortype(ndes,2)=rg
	  colortype(ndes,3)=rb
	  radtype(ndes)=rad
	goto 7020
 7030	write(6,*) ' atom descriptors: ',ndes
	do i=1,ndes
	write(6,*) "type, color, radius ",i,(colortype(i,j),j=1,3),
     &     radtype(i)
	enddo
 21	format(a6,a10)
c --- read atoms and classify ---
	n=0
	nsu=0
	chain=" "
	nbiomat=0
	nbiochain=0
 7040	read(1,7100,end=7009) instring
	if (instring(1:5).eq."MODEL") nsu=nsu+1

	if (instring(12:25).eq."BIOMOLECULE: 1") then
 8010	  read(1,7100) instring
	  if (instring(12:16).eq."APPLY") then
	    ich=43
 8020        if (instring(ich:ich).ne." ") then
	       nbiochain=nbiochain+1
	       biochain(nbiochain)=instring(ich:ich)
	       ich=ich+3
	       goto 8020
	     else
	       continue
	     endif
	    write(6,*) "Chains in Biological Assembly"
	    do i=1,nbiochain
	      write(6,*) i,biochain(i)
	    enddo
	    write(6,*)
	   endif
	  if (instring(14:19).eq."BIOMT1") then
	    nbiomat=nbiomat+1
	    read(instring(20:80),*) ib,(biomat(1,j,nbiomat),j=1,4)
 8030	    format(23x,5f10.6,f15.5)
	  endif
	  if (instring(14:19).eq."BIOMT2") then
c    read(instring,8030) (biomat(2,j,nbiomat),j=1,4)
	    read(instring(20:80),*) ib,(biomat(2,j,nbiomat),j=1,4)
	  endif
	  if (instring(14:19).eq."BIOMT3") then
c    read(instring,8030) (biomat(3,j,nbiomat),j=1,4)
	    read(instring(20:80),*) ib,(biomat(3,j,nbiomat),j=1,4)
	  endif
	  if (instring(14:19).eq."      ") then
	    write(6,*) "Number of BIOMT ",nbiomat
	    do ibio=1,nbiomat
	    do im=1,3
	     write(6,*) "BIOMAT ",ibio,(biomat(im,j,ibio),j=1,4)
	    enddo
	    enddo
	    goto 8099
	  endif
	  goto 8010

	endif
c --- end of BIOMAT
 8099	continue

	if ((instring(1:4).ne.'ATOM').and.
     &      (instring(1:6).ne.'HETATM')) goto 7040

     	read(instring,200) ires
	do ides=1,ndes

 	if (instring(1:6).ne.atomdescriptor(ides)(1:6)) goto 7050

 	 do ia=1,10
 	  if (descriptor(ides)(ia:ia).eq.'-') goto 7060
 	  if (instring(12+ia:12+ia).ne.descriptor(ides)(ia:ia)) goto 7050
 7060 	  continue
	enddo

 	 if ((ires.lt.resrange(1,ides)).or.
     &       (ires.gt.resrange(2,ides))) goto 7050

	 if (radtype(ides).eq.0) goto 7040

c --- check if chain is in biological assembly
	if (nbiochain.ne.0) then
	ibioflag=0
	do ibio=1,nbiochain
	if (instring(22:22).eq.biochain(ibio)) ibioflag=1
	enddo
	if (ibioflag.eq.0) goto 7040
	endif

c --- found an atom to save ---
 	 n=n+1
	 read(instring,300) (coord(i,n),i=1,3)
	 type(n)=ides
C --- assign subunits automatically ---
	chain=instring(22:22)
	if (chain.ne.chainlast) then
	  nsu=nsu+1
	  chainlast=chain
	endif
	su(n)=nsu
	 res(n)=ires

C --- apply biomat ---
cif (nbiomat.gt.0) then
crx=coord(1,n)
cry=coord(2,n)
crz=coord(3,n)
cdo imat=1,nbiomat
csu(n+imat)=su(n)
cres(n+imat)=res(n)
ctype(n+imat)=ides
ccoord(1,n+imat)=rx*biomat(1,1,imat)+
c    &                  ry*biomat(1,2,imat)+
c    &                  rz*biomat(1,3,imat)+
c    &                     biomat(1,4,imat)
ccoord(2,n+imat)=rx*biomat(2,1,imat)+
c    &                  ry*biomat(2,2,imat)+
c    &                  rz*biomat(2,3,imat)+
c    &                     biomat(2,4,imat)
ccoord(3,n+imat)=rx*biomat(3,1,imat)+
c    &                  ry*biomat(3,2,imat)+
c    &                  rz*biomat(3,3,imat)+
c    &                     biomat(3,4,imat)
cenddo
cn=n+nbiomat
cendif

	 goto 7040
 7050 	continue
	enddo
	goto 7040
c --- done reading atoms ---
 7009	write(6,*)' atoms read: ', n, ' from: ',inputfile
	write(6,*) " number of subunits: ",nsu
	write(6,*)' '
 7100	format(a80)
 200	format(22x,i4)
 300	format(30x,3f8.3)
	goto 10
c--------------------------------------------------------------------
c       TRANSLATION
 2	continue
	read(5,*) xtrani,ytrani,ztrani
	xtran=xtran+xtrani
	ytran=ytran+ytrani
	ztran=ztran+ztrani
 	write(6,*) 'translation vector : ',xtran,ytran,ztran
	write(6,*)
	goto 10
c--------------------------------------------------------------------
c	Z ROTATION
 5	call clearmatrix(matrixin)
	read(5,*) angle
	write(6,*) 'z rotation : ',angle
c --minus sign is because original rotations were left-handed! (warning: postdoc code)
	angle=-angle*3.141592/180.
	matrixin(1,1)=cos(angle)
	matrixin(1,2)=-sin(angle)
	matrixin(2,1)=sin(angle)
	matrixin(2,2)=cos(angle)
	call catenate(rm,matrixin)
	goto 10
c--------------------------------------------------------------------
c	Y ROTATION
 4	call clearmatrix(matrixin)
	read(5,*) angle
	write(6,*) 'y rotation : ',angle
	angle=-angle*3.141592/180.
	matrixin(1,1)=cos(angle)
	matrixin(1,3)=sin(angle)
	matrixin(3,1)=-sin(angle)
	matrixin(3,3)=cos(angle)
	call catenate(rm,matrixin)
	goto 10
c--------------------------------------------------------------------
c	X ROTATION
 3	call clearmatrix(matrixin)
	read(5,*) angle
	write(6,*) 'x rotation : ',angle
	angle=-angle*3.141592/180.
	matrixin(2,2)=cos(angle)
	matrixin(2,3)=-sin(angle)
	matrixin(3,2)=sin(angle)
	matrixin(3,3)=cos(angle)
	call catenate(rm,matrixin)
	goto 10
c--------------------------------------------------------------------
c       SCALE
 6	continue
	read(5,*) rscalei
	rscale=rscale*rscalei
	write(6,*) 'scale factor : ',rscale
	write(6,*)
	goto 10
c--------------------------------------------------------------------
c	CENTERING of coordinates
c       three options: CEN(TER) will center on max/min coordinates in x,y,z
c                      AUT(O-CENTER) will center the rotated coordinates
c                        in the frame, and place the uppermost atom
c                        just below the image plane
c		         (done when image processing begins)
c       ROTATIONs are applied before the centering
c       TRANSLATIONs are applied after the centering
 7	continue
 	read(5,101) command
 	autocenter=0
 	if (command.eq.'aut') then
	   autocenter=1
	   goto 10
	endif

	if (command.eq.'cen') then
	   autocenter=2
	   goto 10
	endif

	goto 10
c--------------------------------------------------------------------
c	WORLD parameters that describe the environment
 8	continue

c param: rback(3) -- rgb of background, 0.-1.
c param: rfog(3) -- rgb of fog, 0.-1.
c param: pfogh,pfogl -- fractional intensity of fog color at front and back of molecule
	read (5,*) (rback(i),i=1,3),(rfog(i),i=1,3),pfogh,pfogl
	if (pfogl.lt.0.) then
	   pfogl=-pfogl
	endif

	do i=1,3
	if (rback(i).gt.1.) rback(i)=1.
	if (rfog(i).gt.1.) rfog(i)=1.
        colortype(0,i)=rback(i)
	enddo

	if (pfogh.gt.1.) pfogh=1.
	if (pfogl.gt.1.) pfogl=1.

	write (6,202) ' background inten. :',(rback(i),i=1,3)
	write (6,202) ' fog intensity :    ',(rfog(i),i=1,3)
	write (6,202) ' upper fog percent :',pfogh*100.
	write (6,202) ' lower fog percent :',pfogl*100.
	pfogdiff=pfogh-pfogl

c param: icone -- flag, 0=no soft shadow, other=soft shadow
c param: pcone -- fractional shadowing contribution from each atom (0.0023)
c                 larger=darker shadow
c param: coneangle -- angle of shadowing around each atom (2.0)
c                     larger values give tighter shadowed region
c param: rcone -- shadowing applied only if z values of atoms are
c                 greater than rcone (removes many small shadows in
c                 creases) (1.0 A)
c param: pshadowmax -- maximal shadowing amount (0.7) smaller=darker shadow
	read (5,*) icone,pcone,coneangle,rcone,pshadowmax
	if (icone.ne.0) write (6,202) ' draw conical shadows'

c param: ixsize,iysize -- vertical and horizonal size (pixels)
c        negative values are autosized, padded by the value
	read(5,*) ixsize,iysize
	if (ixsize.gt.3000) ixsize=3000
	if (iysize.gt.3000) iysize=3000
	write(6,*) 'input value for image size', ixsize,iysize
	goto 10
 202	format(1x,a20,1x,8f7.2)
 207   	format(/1x,a20/)
c--------------------------------------------------------------------
c	ILLUSTRATION parameters for outlines
 12	continue
	illustrationflag=1

c param: l_low, l_high -- thresholds for drawing contour outlines (3.,8.)
c	                  higher values give fewer outlines
c                         narrower range gives jaggier outlines
c	                  values are based on the number of pixels in the kernel
c param: ikernel -- kernel for doing derivative, values=1,2,3,4 (4)
c param: l_diff_min, l_diff_max -- range of difference in z-values used for calculation (0.,5. A)
c	                  0-1 gives outlines around every atom, 0-1000 outlines the whole molecule
c param: r_low, r_high -- thresholds for drawing subunit outlines (3.,8.)
c param: g_low, g_high -- thresholds for drawing residue outlines (3.,8.)
c param: resdiff -- residue outlines drawn for atoms with this difference in residue numbers or greater

	read(5,*) l_low,l_high,ikernel,l_diff_min,l_diff_max
	read(5,*) r_low,r_high
	read(5,*) g_low,g_high,resdiff
	write(6,*) 'illustration parameters'
	write(6,*) 'l parameters: ',l_low,l_high
	write(6,*) 'g parameters: ',g_low,g_high
	goto 10
c--------------------------------------------------------------------
 111   	write (6,207) ' *begin calculation*'
c --- if no BIOMT in file, use biomat 1 == identity matrix
	nbiomat=max(nbiomat,1)
c ***** Populate conical shadow table ****
	conemax=50.
	do i=-51,51
	do j=-51,51
	 rtable(i,j)=sqrt(float(i)**2+float(j)**2)
	 if (rtable(i,j).gt.conemax) rtable(i,j)=10000.
	enddo
	enddo
	 rtable(0,0)=10000.
c ***** SCALE RADII *****
	do i=1,ndes
	 radtype(i)=radtype(i)*rscale
	 if (radtype(i).gt.radius_max) radius_max=radtype(i)
 119	enddo
c ***** APPLY AUTOCENTERING and AUTOSIZING, if switched on *****
	if (autocenter.gt.0) then

	xmin=10000.
	xmax=-10000.
 	ymin=10000.
	ymax=-10000.
 	zmin=10000.
	zmax=-10000.

	do ia=1,n
	do ibio=1,nbiomat
c--apply biomat
	rx=coord(1,ia)*biomat(1,1,ibio)+coord(2,ia)*biomat(1,2,ibio)+
     &	   coord(3,ia)*biomat(1,3,ibio)+biomat(1,4,ibio)
	ry=coord(1,ia)*biomat(2,1,ibio)+coord(2,ia)*biomat(2,2,ibio)+
     &	   coord(3,ia)*biomat(2,3,ibio)+biomat(2,4,ibio)
	rz=coord(1,ia)*biomat(3,1,ibio)+coord(2,ia)*biomat(3,2,ibio)+
     &	   coord(3,ia)*biomat(3,3,ibio)+biomat(3,4,ibio)

c--apply rotation matrix
	rx2=rx*rm(1,1)+ry*rm(2,1)+rz*rm(3,1)
	ry2=rx*rm(1,2)+ry*rm(2,2)+rz*rm(3,2)
	rz2=rx*rm(1,3)+ry*rm(2,3)+rz*rm(3,3)
	 xmin=MIN(xmin,rx2)
	 xmax=MAX(xmax,rx2)
	 ymin=MIN(ymin,ry2)
	 ymax=MAX(ymax,ry2)
	 zmin=MIN(zmin,rz2)
	 zmax=MAX(zmax,rz2)

	enddo
	enddo

	write(6,*) 'min coordinates : ',xmin,ymin,zmin
	write(6,*) 'max coordinates : ',xmax,ymax,zmax
	xtranc=-xmin-(xmax-xmin)/2.
	ytranc=-ymin-(ymax-ymin)/2.
	if (autocenter.eq.1) then
 	   ztranc=-zmax-radius_max-1.
	   write(6,*) 'automating centering'
	endif
	if (autocenter.eq.2) then
	   ztranc=-zmin-(zmax-zmin)/2.
	   write(6,*) 'x,y,z centering'
	endif

	write(6,*) 'centering vector : ',xtranc,ytranc,ztranc

	if ((ixsize.le.0).or.(iysize.le.0)) then
	  write(6,*)
	  write(6,*) 'applying autosizing'
	  write(6,*) 'x and y frame width: ',-ixsize,-iysize
	  ixsize=-2.*ixsize+2.*radius_max+(xmax-xmin)*rscale
	  iysize=-2.*iysize+2.*radius_max+(ymax-ymin)*rscale
	endif
 	endif

	ixsize=min(ixsize,3000)
	iysize=min(iysize,3000)

c ***** OPEN OUTPUT FILES *****
	ixsize=int(ixsize/2)*2
	iysize=int(iysize/2)*2
	  write(6,*) 'xsize and ysize: ',ixsize,iysize
	  write(6,*)
	 read(5,113) filename
	 write(6,*) "output pnm filename: ",filename
	 open(8,file=filename,form='formatted')
	 write(8,1003) "P3"
 1003	 format(a2)
	 write(8,1004) iysize,ixsize
	 write(8,1004) 255
C this will write an file with opacities
c open(9,file="opacity.pnm",form='formatted')
c write(9,1003) "P3"
c write(9,1004) iysize,ixsize
c write(9,1004) 255
 1004	format(2i5)
 113	format(a80)
c ***** MAP SPHERICAL SURFACES OVER ATOMS *****
	do ix=1,ixsize
	do iy=1,iysize
 	  pix(ix,iy,1)=0.
 	  pix(ix,iy,2)=0.
 	  pix(ix,iy,3)=0.
 	  pix(ix,iy,4)=0.
	  atom(ix,iy)=0
	  bio(ix,iy)=1
 	  zpix(ix,iy)=-10000.
	enddo
	enddo

	if (n.gt.0) then

c ----- create the spherical shading map for atom types -----
	do irad=1,ndes
	ic=1
	irlim=int(radtype(irad))
	if (irlim.gt.100) then
	  write(6,*) 'atoms radius * scale > 100'
	  stop
	endif

	do ix=-irlim-1,irlim+1
	do iy=-irlim-1,irlim+1
	x=float(ix)
	y=float(iy)
	d=sqrt(x*x+y*y)
	if (d.gt.radtype(irad)) goto 350
	z=sqrt(radtype(irad)**2-d*d)
	sphdat(ic,1)=x
	sphdat(ic,2)=y
	sphdat(ic,3)=z
	ic=ic+1
 350	continue
	enddo
 352	enddo
	numpix=ic-1
c ----- then map spherical surface over atoms of the proper type ------
	icount=0

	do ia=1,n

	if (type(ia).ne.irad) goto 500
	icount=icount+1
	do ibio=1,nbiomat
c--apply biomat
	rx=coord(1,ia)*biomat(1,1,ibio)+coord(2,ia)*biomat(1,2,ibio)+
     &	   coord(3,ia)*biomat(1,3,ibio)+biomat(1,4,ibio)
	ry=coord(1,ia)*biomat(2,1,ibio)+coord(2,ia)*biomat(2,2,ibio)+
     &	   coord(3,ia)*biomat(2,3,ibio)+biomat(2,4,ibio)
	rz=coord(1,ia)*biomat(3,1,ibio)+coord(2,ia)*biomat(3,2,ibio)+
     &	   coord(3,ia)*biomat(3,3,ibio)+biomat(3,4,ibio)

c--apply rotation matrix
	rx2=rx*rm(1,1)+ry*rm(2,1)+rz*rm(3,1)
	ry2=rx*rm(1,2)+ry*rm(2,2)+rz*rm(3,2)
	rz2=rx*rm(1,3)+ry*rm(2,3)+rz*rm(3,3)

c--apply centering vector
	rx2=rx2+xtranc
	ry2=ry2+ytranc
	rz2=rz2+ztranc

c--apply translation and scaling
	rx2=(rx2+xtran)*rscale
	ry2=(ry2+ytran)*rscale
	rz2=(rz2+ztran)*rscale

	if (rz2.lt.0.) then

	do ipix=1,numpix
	x=sphdat(ipix,1)+rx2+float(ixsize)/2.
	y=sphdat(ipix,2)+ry2+float(iysize)/2.
	ix=int(x)
	iy=int(y)
	if ((x.gt.float(ixsize)).or.(x.lt.1.).or.
     &      (y.gt.float(iysize)).or.(y.lt.1.)) goto 510
	z=sphdat(ipix,3)+rz2
	if (z.gt.zpix(ix,iy)) then
	 zpix(ix,iy)=z
	 atom(ix,iy)=ia
	 bio(ix,iy)=ibio
	endif
 510	continue
	enddo

	endif

c -- ibio loop
	enddo

 500	continue

c -- ia loop
	enddo

	write(6,*) icount,' spheres added of type: ',irad

c -- irad loop
 340	enddo

	write(6,*) 'shading maps written into depth buffer'

 	endif

c***** CALCULATE SECOND DERIVATIVE OUTLINES ******
c ---- find maximum and minimum z levels ---
c	(note: I use a value of zpix=-10000. to distinguish background)
	zpix_max=-100000.
	zpix_min=100000.
	do ix=1,ixsize
	do iy=1,iysize
	if (zpix(ix,iy).gt.zpix_max) zpix_max=zpix(ix,iy)
	if ((zpix(ix,iy).ne.-10000.).and.
     &	    (zpix(ix,iy).lt.zpix_min)) zpix_min=zpix(ix,iy)
	enddo
	enddo
	zpix_max=min(zpix_max,0.)
	zpix_spread=zpix_max-zpix_min
	write(6,*) 'zpix_min,zpix_max ',zpix_min,zpix_max

c ***** PROCESSING OF THE IMAGE BEGINS HERE*****
	l_diff_min=l_diff_min*rscale
	l_diff_max=l_diff_max*rscale
	write(6,*) ' Pixel processing beginning '
	do ix=1,ixsize
	do iy=1,iysize
	zpix(ix,iy)=min(zpix(ix,iy),0.)
 2009	enddo
 1009	enddo

	do ix=1,ixsize
	do iy=1,iysize

c ***** CONICAL SHADOW TESTING *****
	pconetot=1.
	if ((icone.ne.0).and.(atom(ix,iy).ne.0)) then
	do i=-50,50,5
	do j=-50,50,5
	  if ((ix+i.gt.0).and.(ix+i.lt.ixsize).and.
     &        (iy+j.gt.0).and.(iy+j.lt.iysize)) then
	  rzdiff=zpix(ix+i,iy+j)-zpix(ix,iy)
	   if (rzdiff.gt.rcone) then
	    if (rtable(i,j)*coneangle.lt.rzdiff+rcone) then
	     pconetot=pconetot-pcone
	    endif
	   endif
	  endif
	enddo
	enddo
	pconetot=max(pconetot,pshadowmax)
	endif
c ***** CALCULATE THE FOG PERCENTAGE (DEPTH CUEING) *****
	pfh=pfogh-(zpix_max-zpix(ix,iy))/zpix_spread*pfogdiff
    	if (zpix(ix,iy).lt.zpix_min) pfh=1.
c ***** CALCULATE OUTLINES *****
	g_opacity=0.
	l_opacity=0.
	if (illustrationflag.ne.0) then
c ***** SUBUNIT OUTLINES *****
c ---- (this replaces original calculation of first derivatives)

	if ((ix.gt.1).and.(ix.lt.ixsize).and.
     &      (iy.gt.1).and.(iy.lt.iysize)) then
c --- calculate subunit 'derivatives'  ---
	g=0.
	r=0.
	do i=-2,2
	do j=-2,2
	 if (abs(i*j).ne.4) then
	  if ((su(atom(ix,iy)).ne.su(atom(ix+i,iy+j))).or.
     &        (bio(ix,iy).ne.bio(ix+i,iy+j))) r=r+1.
          if (abs(res(atom(ix,iy))-
     &              res(atom(ix+i,iy+j))).gt.resdiff) g=g+1.
	 endif
	enddo
	enddo
c --- opacities are 1 for completely opaque ---
	g_opacity=min((g-g_low)/(g_high-g_low),1.)
	r_opacity=min((r-r_low)/(r_high-r_low),1.)
	g_opacity=max(g_opacity,r_opacity,0.)
	endif
c ***** SECOND DERIVATIVE OUTLINES *****
	if ((ix.gt.2).and.(ix.lt.ixsize-1).and.
     &      (iy.gt.2).and.(iy.lt.iysize-1)) then
	rl=0.
	l_opacity_ave=0.
	do ixl=-1,1
	do iyl=-1,1
	ixc=ix+ixl
	iyc=iy+iyl

	if (ikernel.eq.1) then
 	 l(ixl,iyl)=abs(1./3. * (
     &	 -0.8*zpix(ixc-1,iyc-1)-1.*zpix(ixc-1,iyc)-0.8*zpix(ixc-1,iyc+1)-
     &	 1.0* zpix(ixc,iyc-1)+7.2*zpix(ixc,iyc)-1.0*zpix(ixc,iyc+1)-
     &	 0.8* zpix(ixc+1,iyc-1)-1.*zpix(ixc+1,iyc)-0.8*zpix(ixc+1,iyc+1)
     &     ))
	endif

	if (ikernel.eq.2) then
	 l(ixl,iyl)=abs(1./3. * (
     &	-0.8*zpix(ixc-1,iyc-1)-1.0*zpix(ixc-1,iyc)-0.8*zpix(ixc-1,iyc+1)-
     &	 1.0*zpix(ixc,iyc-1)+8.8*zpix(ixc,iyc)-1.0*zpix(ixc,iyc+1)-
     &	 0.8*zpix(ixc+1,iyc-1)-1.0*zpix(ixc+1,iyc)-0.8*zpix(ixc+1,iyc+1)-
     &	 0.1*zpix(ixc+2,iyc-1)-0.2*zpix(ixc+2,iyc)-0.1*zpix(ixc+2,iyc+1)-
     &	 0.1*zpix(ixc-2,iyc-1)-0.2*zpix(ixc-2,iyc)-0.1*zpix(ixc-2,iyc+1)-
     &	 0.1*zpix(ixc-1,iyc+2)-0.2*zpix(ixc,iyc+2)-0.1*zpix(ixc+1,iyc+2)-
     &	 0.1*zpix(ixc-1,iyc-2)-0.2*zpix(ixc,iyc-2)-0.1*zpix(ixc+1,iyc-2)
     &     ))
	endif

	if (ikernel.eq.3) then
	do i=-1,1
	do j=-1,1
		rd=abs(zpix(ix,iy)-zpix(ix+i,iy+j))
		if (rd.gt.l_diff_min) then
		rd=(rd-l_diff_min)/(l_diff_max-l_diff_min)
		l(ixl,iyl)=l(ixl,iyl)+min(rd,1.)
		endif
	enddo
	enddo
	endif

	if (ikernel.eq.4) then
	do i=-2,2
	do j=-2,2
	 if (abs(i*j).ne.4) then
		rd=abs(zpix(ix,iy)-zpix(ix+i,iy+j))
		if (rd.gt.l_diff_min) then
		rd=abs(zpix(ix,iy)-zpix(ix+i,iy+j))
		rd=(rd-l_diff_min)/(l_diff_max-l_diff_min)
		l(ixl,iyl)=l(ixl,iyl)+min(rd,1.)
		endif
	 endif
	enddo
	enddo
	endif

	l(ixl,iyl)=min((l(ixl,iyl)-l_low)/(l_high-l_low),1.)
	l(ixl,iyl)=max(l(ixl,iyl),0.)
	if (l(ixl,iyl).gt.0.) rl=rl+1.
	l_opacity_ave=l_opacity_ave+l(ixl,iyl)
	enddo
	enddo
	if (rl.ge.6.) then
	l_opacity=l_opacity_ave/6.
	else
	l_opacity=l(0,0)
	endif
	l_opacity=min(l_opacity,1.)
	l_opacity=max(l_opacity,0.)
	endif
c ---- combine subunit outlines and derivative outlines ----
	l_opacity=max(l_opacity,g_opacity)
	endif
c ***** CALCULATE THE TOTAL PIXEL INTENSITY *****
	ropacity=0.
	do icolor=1,3

	rcolor=
     &     pfh*(pconetot*(colortype(type(atom(ix,iy)),icolor)))+
     &     (1.-pfh)*rfog(icolor)

	pix(ix,iy,icolor)=(1.-l_opacity)*rcolor

c ----calculate pixel opacity
	if (type(atom(ix,iy)).ne.0) ropacity=1.
	pix(ix,iy,4)=max(ropacity,l_opacity)

	enddo

 	enddo

c ***** output of a scan line *****
c ----- PPM format -----
	iscan=0
	do iout=1,iysize
	do ic=1,3
	iscan=iscan+1
	scanline(iscan)=int(pix(ix,iout,ic)*255.)
	scanline(iscan)=min(scanline(iscan),255)
	scanline(iscan)=max(scanline(iscan),0)
	enddo
	enddo
	write(8,1002) (scanline(if),if=1,iysize*3)
C -- write opacity
ciscan=0
cdo iout=1,iysize
cdo ic=1,3
ciscan=iscan+1
cscanline(iscan)=int(pix(ix,iout,4)*255.)
cscanline(iscan)=min(scanline(iscan),255)
cscanline(iscan)=max(scanline(iscan),0)
cenddo
cenddo
cwrite(9,1002) (scanline(if),if=1,iysize*3)
 1002	format(20i4)
c ----- diagnostic ------
	if (int(ix/20)*20.eq.int((float(ix)/20.)*20.)) then
	 write(6,6669) (int(pix(ix,iyo,1)*9.),iyo=1,iysize,20)
 6669	format(65i1)
	endif

 1000	enddo

 999	stop
	end
c--------------------------------------------------------------------
	subroutine catenate(m1,m2)
	real*4 m1(4,4),m2(4,4),m(4,4)
c--------------------------------------------------------------------
	m(1,1)=m1(1,1)*m2(1,1)+m1(2,1)*m2(1,2)+m1(3,1)*m2(1,3)+
     &	       m1(4,1)*m2(1,4)
	m(1,2)=m1(1,2)*m2(1,1)+m1(2,2)*m2(1,2)+m1(3,2)*m2(1,3)+
     &	       m1(4,2)*m2(1,4)
	m(1,3)=m1(1,3)*m2(1,1)+m1(2,3)*m2(1,2)+m1(3,3)*m2(1,3)+
     &	       m1(4,3)*m2(1,4)
	m(1,4)=m1(1,4)*m2(1,1)+m1(2,4)*m2(1,2)+m1(3,4)*m2(1,3)+
     &	       m1(4,4)*m2(1,4)
	m(2,1)=m1(1,1)*m2(2,1)+m1(2,1)*m2(2,2)+m1(3,1)*m2(2,3)+
     &	       m1(4,1)*m2(2,4)
	m(2,2)=m1(1,2)*m2(2,1)+m1(2,2)*m2(2,2)+m1(3,2)*m2(2,3)+
     &	       m1(4,2)*m2(2,4)
	m(2,3)=m1(1,3)*m2(2,1)+m1(2,3)*m2(2,2)+m1(3,3)*m2(2,3)+
     &	       m1(4,3)*m2(2,4)
	m(2,4)=m1(1,4)*m2(2,1)+m1(2,4)*m2(2,2)+m1(3,4)*m2(2,3)+
     &	       m1(4,4)*m2(2,4)
	m(3,1)=m1(1,1)*m2(3,1)+m1(2,1)*m2(3,2)+m1(3,1)*m2(3,3)+
     &	       m1(4,1)*m2(3,4)
	m(3,2)=m1(1,2)*m2(3,1)+m1(2,2)*m2(3,2)+m1(3,2)*m2(3,3)+
     &	       m1(4,2)*m2(3,4)
	m(3,3)=m1(1,3)*m2(3,1)+m1(2,3)*m2(3,2)+m1(3,3)*m2(3,3)+
     &	       m1(4,3)*m2(3,4)
	m(3,4)=m1(1,4)*m2(3,1)+m1(2,4)*m2(3,2)+m1(3,4)*m2(3,3)+
     &	       m1(4,4)*m2(3,4)
	m(4,1)=m1(1,1)*m2(4,1)+m1(2,1)*m2(4,2)+m1(3,1)*m2(4,3)+
     &	       m1(4,1)*m2(4,4)
	m(4,2)=m1(1,2)*m2(4,1)+m1(2,2)*m2(4,2)+m1(3,2)*m2(4,3)+
     &	       m1(4,2)*m2(4,4)
	m(4,3)=m1(1,3)*m2(4,1)+m1(2,3)*m2(4,2)+m1(3,3)*m2(4,3)+
     &	       m1(4,3)*m2(4,4)
	m(4,4)=m1(1,4)*m2(4,1)+m1(2,4)*m2(4,2)+m1(3,4)*m2(4,3)+
     &	       m1(4,4)*m2(4,4)
	do j=1,4
	do i=1,4
 	m1(i,j)=m(i,j)
	enddo
	enddo
	return
	end
c--------------------------------------------------------------------
	subroutine clearmatrix(m)
	real*4 m(4,4)
c--------------------------------------------------------------------
	do i=1,4
	do j=1,4
	x=0.
	if (i.eq.j) x=1.
1	m(i,j)=x
	enddo
	enddo
	return
	end

max_issues_repo_name: shozebhaider/Illustrate

max_forks_repo_name: shozebhaider/Illustrate


text:
// Generated Files
${PROJ_DIR}/lpif/rtl/lpif_txrx/lpif_txrx_x8_f2_master_top.sv
${PROJ_DIR}/lpif/rtl/lpif_txrx/lpif_txrx_x8_f2_master_concat.sv
${PROJ_DIR}/lpif/rtl/lpif_txrx/lpif_txrx_x8_f2_master_name.sv

// Logic Link files
-f ${PROJ_DIR}/llink/rtl/llink.f

// Common Files
-f ${PROJ_DIR}/common/rtl/common.f

ext: f

max_stars_repo_name: chipsalliance/aib-protocols

max_issues_repo_name: chipsalliance/aib-protocols

max_forks_repo_name: chipsalliance/aib-protocols


text:
subroutine euler2a(neqn,t0,tf,u0,nsteps,u)
!
! Subroutine euler2a computes an ODE solution by a fixed
! step modified Euler method for a series of points along
! the solution by repeatedly calling subroutine sseuler for
! a single modified Euler step.
!
! Argument list
!
! neqn number of first order ODEs
!
! t0 initial value of independent variable
!
! tf final value of independent variable
!
! u0 initial condition vector of length neqn
!
! nsteps number of modified Euler steps
!
! u ODE solution vector of length neqn after
! nsteps steps
!
! Double precision coding is used
implicit double precision(a-h,o-z)
!
! Size arrays
dimension u0(neqn), u(neqn), e(neqn)
!
! Integration step
h=(tf-t0)/dfloat(nsteps)
!
! nsteps modified Euler steps
do j=1,nsteps
!
! Modified Euler step
call sseuler(neqn,t0,u0,h,t,u,e)
!
! Reset base point values for the next modified
! Euler step
do i=1,neqn
u0(i)=u(i)
end do
t0=t
!
! Next modified Euler step
end do
!
! nsteps modified Euler steps completed
return
!
! End of euler2a
end

max_stars_repo_name: cunyizju/Runge-Kutta-ODE-FORTRAN

max_issues_repo_name: cunyizju/Runge-Kutta-ODE-FORTRAN

max_forks_repo_name: cunyizju/Runge-Kutta-ODE-FORTRAN


    program demo_btest
    implicit none
    integer :: i, j, pos, a(2,2)
    logical :: bool
    character(len=*),parameter :: g='(*(g0))'

         i = 32768 + 1024 + 64
        write(*,'(a,i0,"=>",b32.32,/)')'Looking at the integer: ',i

        ! looking one bit at a time from LOW BIT TO HIGH BIT
        write(*,g)'from bit 0 to bit ',bit_size(i),'==>'
        do pos=0,bit_size(i)-1
            bool = btest(i, pos)
            write(*,'(l1)',advance='no')bool
        enddo
        write(*,*)

        ! a binary format the hard way.
        ! Note going from bit_size(i) to zero.
        write(*,*)
        write(*,g)'so for ',i,' with a bit size of ',bit_size(i)
        write(*,'(b32.32)')i
        write(*,g)merge('^','_',[(btest(i,j),j=bit_size(i)-1,0,-1)])
        write(*,*)
        write(*,g)'and for ',-i,' with a bit size of ',bit_size(i)
        write(*,'(b32.32)')-i
        write(*,g)merge('^','_',[(btest(-i,j),j=bit_size(i)-1,0,-1)])

        ! elemental:
        !
        a(1,:)=[ 1, 2 ]
        a(2,:)=[ 3, 4 ]
        write(*,*)
        write(*,'(a,/,*(i2,1x,i2,/))')'given the array a ...',a
        ! the second bit of all the values in a
        write(*,'(a,/,*(l2,1x,l2,/))')'the value of btest (a, 2)',btest(a,2)
        ! bits 1,2,3,4 of the value 2
        write(*,'(a,/,*(l2,1x,l2,/))')'the value of btest (2, a)',btest(2,a)
    end program demo_btest


text:
module mod_chkdiv
  use mpi
  use mod_common_mpi, only: myid,coord,ierr
  use mod_types
  implicit none
  private
  public chkdiv
  contains
  subroutine chkdiv(n,dli,dzfi,u,v,w,divtot,divmax)
    !
    ! checks the divergence of the velocity field
    !
    implicit none
    integer , intent(in), dimension(3) :: n
    real(rp), intent(in), dimension(3) :: dli
    real(rp), intent(in), dimension(0:) :: dzfi
    real(rp), intent(in), dimension(0:,0:,0:) :: u,v,w
    real(rp), intent(out) :: divtot,divmax
    real(rp) :: dxi,dyi,div!,dzi,div
    integer :: i,j,k,im,jm,km
    !integer :: ii,jj
    !
    dxi = dli(1)
    dyi = dli(2)
    !dzi = dli(3)
    divtot = 0.
    divmax = 0.
    !$OMP PARALLEL DO DEFAULT(none) &
    !$OMP SHARED(n,u,v,w,dxi,dyi,dzfi) &
    !$OMP PRIVATE(i,j,k,im,jm,km,div) &
    !$OMP REDUCTION(+:divtot) &
    !$OMP REDUCTION(max:divmax)
    do k=1,n(3)
       km = k-1
       do j=1,n(2)
          jm = j-1
          do i=1,n(1)
             im = i-1
             div = (w(i,j,k)-w(i,j,km))*dzfi(k) + &
                   (v(i,j,k)-v(i,jm,k))*dyi     + &
                   (u(i,j,k)-u(im,j,k))*dxi
             divmax = max(divmax,abs(div))
             divtot = divtot + div
             !ii = coord(1)*n(1)+i
             !jj = coord(2)*n(2)+j
             !if(abs(div).ge.1.e-12) print*,div,'Large divergence at grid cell: ',ii,jj,k,div
          enddo
       enddo
    enddo
    !$OMP END PARALLEL DO
    call mpi_allreduce(MPI_IN_PLACE,divtot,1,MPI_REAL_RP,MPI_SUM,MPI_COMM_WORLD,ierr)
    call mpi_allreduce(MPI_IN_PLACE,divmax,1,MPI_REAL_RP,MPI_MAX,MPI_COMM_WORLD,ierr)
    if(myid.eq.0) print*, 'Total divergence = ', divtot, '| Maximum divergence = ', divmax
    return
  end subroutine chkdiv
end module mod_chkdiv

max_stars_repo_name: zhoutengye/VSIAM-fft

max_issues_repo_name: kaolagithub/CaNS

max_forks_repo_name: kaolagithub/CaNS


text: Thanks for adding to the featured page!  Please also help add to the realdeal page for the topic:  Noise Ordinance.   Users/PhilipNeustrom

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!----------------------------------------------------------------
! (c) Copyright, 2016 by the Regents of the University of California.
! DTLclass: Drift-tube-linac beam line element class
!             in Lattice module of APPLICATION layer.
! Version: 1.0
! Author: Ji Qiang, LBNL
! Description: This class defines the linear transfer map and RF field
!              for the DTL beam line elment.
! Comments:
!----------------------------------------------------------------
      module DTLclass
        use PhysConstclass
        use Dataclass
        integer, private, parameter :: Nparam = 25
        type DTL
          !Itype = 101
          integer :: Nseg,Mapstp,Itype
          double precision :: Length
          double precision, dimension(Nparam) :: Param
          ! Param(1) : zedge
          !      (2) : scale
          !      (3) : RF frequency
          !      (4) : theta0
          !      (5) : file ID
          !      (6) : radius
          !      (7) : quad 1 length
          !      (8) : quad 1 gradient
          !      (9) : quad 2 length
          !      (10) : quad 2 gradient
          !      (11) : x misalignment error for Quad 1
          !      (12) : y misalignment error for Quad 1
          !      (13) : rotation error x for Quad 1
          !      (14) : rotation error y for Quad 1
          !      (15) : rotation error z for Quad 1
          !      (16) : x misalignment error for Quad 2
          !      (17) : x misalignment error for Quad 2
          !      (18) : rotation error x for Quad 2
          !      (19) : rotation error y for Quad 2
          !      (20) : rotation error z for Quad 2
          !      (21) : x misalignment error for RF cavity
          !      (22) : y misalignment error for RF cavity
          !      (23) : rotation error x for RF cavity
          !      (24) : rotation error y for RF cavity
          !      (25) : rotation error z for RF cavity
        end type DTL
        interface getparam_DTL
          module procedure getparam1_DTL,  &
                          getparam2_DTL,   &
                          getparam3_DTL
        end interface
        interface setparam_DTL
          module procedure setparam1_DTL,  &
                          setparam2_DTL, setparam3_DTL
        end interface
      contains
        subroutine construct_DTL(this,numseg,nmpstp,type,blength)
        implicit none
        type (DTL), intent(out) :: this
        integer, intent(in) :: numseg,nmpstp,type
        double precision, intent(in) :: blength

        this%Nseg = numseg
        this%Mapstp = nmpstp
        this%Itype = type
        this%Length = blength

        this%Param = 0.0

        end subroutine construct_DTL

        subroutine setparam1_DTL(this,i,value)
        implicit none
        type (DTL), intent(inout) :: this
        integer, intent(in) :: i
        double precision, intent(in) :: value

        this%Param(i) = value

        end subroutine setparam1_DTL

        subroutine setparam2_DTL(this,values)
        implicit none
        type (DTL), intent(inout) :: this
        double precision, dimension(:), intent(in) :: values

        this%Param = values

        end subroutine setparam2_DTL

        subroutine setparam3_DTL(this,numseg,nmpstp,type,blength)
        implicit none
        type (DTL), intent(inout) :: this
        integer, intent(in) :: numseg,nmpstp,type
        double precision, intent(in) :: blength

        this%Nseg = numseg
        this%Mapstp = nmpstp
        this%Itype = type
        this%Length = blength

        end subroutine setparam3_DTL

        subroutine getparam1_DTL(this,i,blparam)
        implicit none
        type (DTL), intent(in) :: this
        integer, intent(in) :: i
        double precision, intent(out) :: blparam

        blparam = this%Param(i)

        end subroutine getparam1_DTL

        subroutine getparam2_DTL(this,blparams)
        implicit none
        type (DTL), intent(in) :: this
        double precision, dimension(:), intent(out) :: blparams

        blparams = this%Param

        end subroutine getparam2_DTL

        subroutine getparam3_DTL(this,blength,bnseg,bmapstp,&
                                       btype)
        implicit none
        type (DTL), intent(in) :: this
        double precision, intent(out) :: blength
        integer, intent(out) :: bnseg,bmapstp,btype

        blength = this%Length
        bnseg = this%Nseg
        bmapstp = this%Mapstp
        btype = this%Itype

        end subroutine getparam3_DTL

        subroutine maplinear_DTL(t,tau,xm,this,refpt,Bchg,Bmass)
        implicit none
        include 'mpif.h'
        double precision, intent(in) :: t,tau,Bchg,Bmass
        double precision, dimension(6,6), intent(out) :: xm
        double precision, dimension(6), intent(inout) :: refpt
        type (DTL), intent(in) :: this
        double precision, dimension(18) :: y
        double precision :: h,t0,gi,betai,ui,squi,sqi3
        double precision :: dlti,thli,gf,betaf,uf,squf,sqf3
        double precision :: dltf,thlf,zedge,escale,ww,theta0,qmcc
        double precision :: ein,e1in,e2in,ef,e1f,e2f
        double precision :: sinthi,costhi,uprimi,qpwi,tfin,sinthf,costhf
        double precision :: uprimf,qpwf
        integer :: mpstp

        xm = 0.0

        zedge = this%Param(1)
        escale = this%Param(2)
        ww = this%Param(3)/Scfreq
        theta0 = this%Param(4)*asin(1.0)/90
        qmcc = Bchg/Bmass
        mpstp = this%Mapstp

        y(1)=0.0
        y(2)=0.0
        y(3)=0.0
        y(4)=0.0
        y(5)=refpt(5)
        y(6)=refpt(6)
        y(7)=1.0
        y(8)=0.0
        y(9)=0.0
        y(10)=1.0
        y(11)=1.0
        y(12)=0.0
        y(13)=0.0
        y(14)=1.0
        y(15)=1.0
        y(16)=0.0
        y(17)=0.0
        y(18)=1.0
        h=tau/mpstp
        t0=t
        gi=-refpt(6)
        betai=sqrt(gi**2-1.)/gi
        ui=gi*betai
        squi=sqrt(ui)
        sqi3=sqrt(ui)**3

        call getaxfldE_DTL(t,this,ein,e1in,e2in)
        sinthi=sin(ww*refpt(5)+theta0)
        costhi=cos(ww*refpt(5)+theta0)
        uprimi=qmcc*ein/betai*costhi
        qpwi=0.5*qmcc*Scxl/(ui*ww)
        dlti=Scxl*(0.5*uprimi/ui-qpwi*e1in*sinthi)
        thli=1.5*Scxl*(uprimi/ui)

        call rk6i_DTL(h,mpstp,t0,y,18,this,Bchg,Bmass)

        refpt(5)=y(5)
        refpt(6)=y(6)
        gf=-refpt(6)
        betaf=sqrt(gf**2-1.)/gf
        uf=gf*betaf
        squf=sqrt(uf)
        sqf3=sqrt(uf)**3
        tfin = t + tau
        call getaxfldE_DTL(tfin,this,ef,e1f,e2f)
        sinthf=sin(ww*refpt(5)+theta0)
        costhf=cos(ww*refpt(5)+theta0)
        uprimf=qmcc*ef/betaf*costhf
        qpwf=0.5*qmcc*Scxl/(uf*ww)
        dltf=Scxl*(0.5*uprimf/uf-qpwf*e1f*sinthf)
        thlf=1.5*Scxl*(uprimf/uf)

        xm(1,1)= y(7)*squi/squf+y(9)*squi/squf*dlti
        xm(2,1)=(y(8)-y(7)*dltf)*squi*squf+(y(10)-y(9)*dltf)*squi*squf*&
                 dlti
        xm(1,2)= y(9)/(squi*squf)
        xm(2,2)=(y(10)-y(9)*dltf)*squf/squi
        xm(3,3)= y(11)*squi/squf+y(13)*squi/squf*dlti
        xm(4,3)= &
        (y(12)-y(11)*dltf)*squi*squf+(y(14)-y(13)*dltf)*squi*squf*dlti
        xm(3,4)= y(13)/(squi*squf)
        xm(4,4)=(y(14)-y(13)*dltf)*squf/squi
        xm(5,5)= y(15)*sqi3/sqf3+y(17)*sqi3/sqf3*thli
        xm(6,5)= &
        (y(16)-y(15)*thlf)*sqi3*sqf3+(y(18)-y(17)*thlf)*sqi3*sqf3*thli
        xm(5,6)= y(17)/(sqi3*sqf3)
        xm(6,6)=(y(18)-y(17)*thlf)*sqf3/sqi3

        !print*,"xm: ",xm(1,1),xm(2,1),xm(1,2),xm(2,2),xm(3,3),xm(4,3),&
        !        xm(4,4),xm(5,5),xm(5,6),xm(6,6)

        end subroutine maplinear_DTL

        subroutine rk6i_DTL(h,ns,t,y,nvar,this,Bchg,Bmass)
        implicit none
        include 'mpif.h'
        integer, intent(in) :: ns,nvar
        double precision, intent(inout) :: t
        double precision, intent(in) :: h,Bchg,Bmass
        double precision, dimension(nvar), intent(inout) :: y
        type (DTL), intent(in) :: this
        double precision, dimension(nvar) :: yt,a,b,c,d,e,f,g,o,p
        double precision:: tint,tt
        integer :: i
        tint=t
        do i=1,ns
          call intfunc1_DTL(t,y,f,this,Bchg,Bmass)
          a=h*f
          yt=y+a/9.d0
          tt=t+h/9.d0
          call intfunc1_DTL(tt,yt,f,this,Bchg,Bmass)
          b=h*f
          yt=y + (a + 3.d0*b)/24.d0
          tt=t+h/6.d0
          call intfunc1_DTL(tt,yt,f,this,Bchg,Bmass)
          c=h*f
          yt=y+(a-3.d0*b+4.d0*c)/6.d0
          tt=t+h/3.d0
          call intfunc1_DTL(tt,yt,f,this,Bchg,Bmass)
          d=h*f
          yt=y + (278.d0*a - 945.d0*b + 840.d0*c + 99.d0*d)/544.d0
          tt=t+.5d0*h
          call intfunc1_DTL(tt,yt,f,this,Bchg,Bmass)
          e=h*f
          yt=y + (-106.d0*a+273.d0*b-104.d0*c-107.d0*d+48.d0*e)/6.d0
          tt = t+2.d0*h/3.d0
          call intfunc1_DTL(tt,yt,f,this,Bchg,Bmass)
          g=h*f
          yt = y+(110974.d0*a-236799.d0*b+68376.d0*c+ &
                  103803.d0*d-10240.d0*e + 1926.d0*g)/45648.d0
          tt = t + 5.d0*h/6.d0
          call intfunc1_DTL(tt,yt,f,this,Bchg,Bmass)
          o=h*f
          yt = y+(-101195.d0*a+222534.d0*b-71988.d0*c-  &
                  26109.d0*d-2.d4*e-72.d0*g+22824.d0*o)/25994.d0
          tt = t + h
          call intfunc1_DTL(tt,yt,f,this,Bchg,Bmass)
          p=h*f
          y = y+(41.d0*a+216.d0*c+27.d0*d+ &
                 272.d0*e+27.d0*g+216.d0*o+41.d0*p)/840.d0
          t=tint+i*h
        enddo

        end subroutine rk6i_DTL

        subroutine intfunc1_DTL(t,y,f,this,Bchg,Bmass)
        implicit none
        include 'mpif.h'
        double precision, intent(in) :: t,Bchg,Bmass
        double precision, dimension(:), intent(in) :: y
        double precision, dimension(:), intent(out) :: f
        type (DTL), intent(in) :: this
        double precision :: gamma0,beta0,gbet,s11,s33,s55,s11tmp
        double precision :: zedge,escale,ez1,ezp1,ezpp1,ww,theta0,qmcc,&
                            sinphi,cosphi,rfdsgn,brho,bgrad
        integer :: my_rank, ierr

        zedge = this%Param(1)
        call getaxfldE_DTL(t,this,ez1,ezp1,ezpp1)
        call getBgradfld_DTL(t,this,bgrad)
        ww = this%Param(3)/Scfreq
        theta0 = this%Param(4)*asin(1.0)/90
        qmcc = Bchg/Bmass

        f(1) = 0.0
        f(2) = 0.0
        f(3) = 0.0
        f(4) = 0.0

        ! synchronous particle:
        gamma0=-y(6)
        beta0=sqrt(gamma0**2-1.)/gamma0
        gbet=sqrt((gamma0-1.0)*(gamma0+1.0))
        sinphi=sin(ww*y(5)+theta0)
        cosphi=cos(ww*y(5)+theta0)
        rfdsgn=ez1*cosphi
        f(5)=1.0/(beta0*Scxl)
        f(6)=-qmcc*rfdsgn

        ! matrix elements
        brho=gbet/Clight/qmcc
        s11tmp=0.5e0*(1.0+0.5*gamma0**2)* &
         (qmcc*rfdsgn/beta0**2/gamma0**2)**2 + &
          qmcc*0.5/beta0**3/gamma0**3*ez1*sinphi*ww/Scxl
        s11=(s11tmp+bgrad/brho)*Scxl
        s33=(s11tmp-bgrad/brho)*Scxl
        s55=  &
         -1.5e0*qmcc/beta0**2/gamma0*ezp1*cosphi &
         +(beta0**2+0.5)/beta0**3/gamma0*qmcc*ez1*sinphi*ww/Scxl &
         +1.5e0*(1.0-0.5*gamma0**2)* &
         (qmcc*rfdsgn/beta0**2/gamma0**2)**2
        s55=Scxl*s55

        f(7)=y(8)/Scxl
        f(8)=-s11*y(7)
        f(9)=y(10)/Scxl
        f(10)=-s11*y(9)
        f(11)=y(12)/Scxl
        f(12)=-s33*y(11)
        f(13)=y(14)/Scxl
        f(14)=-s33*y(13)
        f(15)=y(16)/Scxl
        f(16)=-s55*y(15)
        f(17)=y(18)/Scxl
        f(18)=-s55*y(17)

        end subroutine intfunc1_DTL

        !interpolate the field from the DTL rf cavity onto bunch location.
        subroutine getaxfldE_DTL(z,this,ez1,ezp1,ezpp1)
        implicit none
        include 'mpif.h'
        type (DTL), intent(in) :: this
        double precision, intent(in) :: z
        double precision, intent(out) :: ez1,ezp1,ezpp1
        double precision:: zz,hstep,slope,zedge,escale
        integer :: klo,khi,k
        integer :: my_rank,ierr

        zedge = this%Param(1)
        escale = this%Param(2)
        zz=z-zedge
        klo=1
        khi=Ndata
1       if(khi-klo.gt.1) then
          k=(khi+klo)/2
          if(zdat(k).gt.zz)then
             khi=k
          else
             klo=k
          endif
          goto 1
        endif
        hstep=zdat(khi)-zdat(klo)
        slope=(edat(khi)-edat(klo))/hstep
        ez1 =edat(klo)+slope*(zz-zdat(klo))
        slope=(epdat(khi)-epdat(klo))/hstep
        ezp1=epdat(klo)+slope*(zz-zdat(klo))
        ezpp1 = 0.0
        ez1 = ez1*escale
        ezp1 = ezp1*escale

        !print*,"ez1: ",ez1,escale,zedge,zz
!        call MPI_COMM_RANK(MPI_COMM_WORLD,my_rank,ierr)
!        if(my_rank.eq.1) then
!         write(19,101)z,ez1/1.0e6,ezp1/1.0e6,zz,float(m),float(init)
!        endif
!        if(my_rank.eq.1) then
!         write(19,101)z,gt
!        endif
  101   format(6(1x,1pe13.6))

        end subroutine getaxfldE_DTL

        subroutine getBgradfld_DTL(z,this,bgrad)
        implicit none
        include 'mpif.h'
        type (DTL), intent(in) :: this
        double precision, intent(in) :: z
        double precision, intent(out) :: bgrad
        double precision:: zz,bgrad1,bgrad2,len1,len2,len,zedge

        len = this%Length
        zedge = this%Param(1)
        len1 = this%Param(7)
        bgrad1 = this%Param(8)
        len2 = len - this%Param(9)
        bgrad2 = this%Param(10)
        zz=z-zedge
        if(zz.lt.len1) then
          bgrad = bgrad1
        else if(zz.lt.len2) then
          bgrad = 0.0
        else
          bgrad = bgrad2
        endif

        end subroutine getBgradfld_DTL

        !get external field with displacement and rotation errors.
        subroutine  getflderrold_DTL(pos,extfld,this,dx,dy,anglex,angley,&
                                  anglez)
        implicit none
        include 'mpif.h'
        double precision, dimension(4), intent(in) :: pos
        type (DTL), intent(in) :: this
        double precision, intent(in) :: dx,dy,anglex,angley,anglez
        double precision, dimension(6), intent(out) :: extfld
        double precision :: zedge,escale,theta0,ww,len,tt,zz,xl
        double precision :: ez1,ezp1,ezpp1,ezppp,f1,f1p,clite,tmpsin,tmpcos,pi
        double precision :: tmpex,tmpey,tmpez,tmpbx,tmpby,tmpbz
        double precision :: bgrad,bgrad1,bgrad2,len1,len2,r2,zmid
        double precision :: dx1,dy1,anglex1,angley1,anglez1
        double precision, dimension(3) :: temp,tmp
        integer :: i

        clite = 299792458.e0
        pi = 2*asin(1.0)

        zedge = this%Param(1)
        escale = this%Param(2)
        len = this%Length
        zmid = zedge + len/2
        ww = this%Param(3)/Scfreq
        xl = Scxl/ww  !real frequency has to be used here
        theta0 = this%Param(4)*asin(1.0)/90
        !move into the tank local coordinate.
        !zz=pos(3)-zedge

        len1 = this%Param(7)
        bgrad1 = this%Param(8)
        len2 = len - this%Param(9)
        bgrad2 = this%Param(10)

!        dx = this%Param(11)
!        dy = this%Param(12)
!        anglex = this%Param(13)
!        angley = this%Param(14)
!        anglez = this%Param(15)
        temp(1) = pos(1) - dx
        temp(2) = pos(2) - dy
        tmp(1) = temp(1)*cos(anglez) + temp(2)*sin(anglez)
        tmp(2) = -temp(1)*sin(anglez) + temp(2)*cos(anglez)
        tmp(3) = pos(3) - zedge
        temp(1) = tmp(1)*cos(angley)+tmp(3)*sin(angley)
        temp(2) = tmp(2)
        temp(3) = -tmp(1)*sin(angley)+tmp(3)*cos(angley)
        tmp(1) = temp(1)
        tmp(2) = temp(2)*cos(anglex)+temp(3)*sin(anglex)
        tmp(3) = -temp(2)*sin(anglex)+temp(3)*cos(anglex)
        zz = tmp(3)

        if(zz.lt.len1) then
          bgrad = bgrad1
        else if(zz.lt.len2) then
          bgrad = 0.0
        else
          bgrad = bgrad2
        endif
        extfld(1) = 0.0
        extfld(2) = 0.0
        extfld(3) = 0.0
        extfld(4) = bgrad*tmp(2)
        extfld(5) = bgrad*tmp(1)
        extfld(6) = 0.0

        tmp(1) = extfld(4)
        tmp(2) = extfld(5)*cos(anglex)-extfld(6)*sin(anglex)
        tmp(3) = extfld(5)*sin(anglex)+extfld(6)*cos(anglex)
        temp(1) = tmp(1)*cos(angley)-tmp(3)*sin(angley)
        temp(2) = tmp(2)
        temp(3) = tmp(1)*sin(angley)+tmp(3)*cos(angley)
        extfld(4) = temp(1)*cos(anglez) - temp(2)*sin(anglez)
        extfld(5) = temp(1)*sin(anglez) + temp(2)*cos(anglez)
        extfld(6) = temp(3)

! get external RF field on axis from analytical function
!-----------------------------------------------------------------
        dx1 = this%Param(16)
        dy1 = this%Param(17)
        anglex1 = this%Param(18)
        angley1 = this%Param(19)
        anglez1 = this%Param(20)

        zz=pos(3)-zmid
        temp(1) = pos(1) - dx
        temp(2) = pos(2) - dy
        tmp(1) = temp(1)*cos(anglez) + temp(2)*sin(anglez)
        tmp(2) = -temp(1)*sin(anglez) + temp(2)*cos(anglez)
        tmp(3) = zz
        temp(1) = tmp(1)*cos(angley)+tmp(3)*sin(angley)
        temp(2) = tmp(2)
        temp(3) = -tmp(1)*sin(angley)+tmp(3)*cos(angley)
        tmp(1) = temp(1)
        tmp(2) = temp(2)*cos(anglex)+temp(3)*sin(anglex)
        tmp(3) = -temp(2)*sin(anglex)+temp(3)*cos(anglex)
        zz = tmp(3)

        ez1 = Fcoef(1)/2
        ezp1 = 0.0
        ezpp1 = 0.0
        ezppp = 0.0
        do i = 2, (Ndata-1)/2 + 1
          ez1 = ez1 + Fcoef(2*i-2)*cos((i-1)*2*pi*zz/len) + &
                Fcoef(2*i-1)*sin((i-1)*2*pi*zz/len)
          ezp1 = ezp1 + (i-1)*2*pi/len*(-Fcoef(2*i-2)*sin((i-1)*2*pi*zz/len)+&
                Fcoef(2*i-1)*cos((i-1)*2*pi*zz/len))
          ezpp1 = ezpp1+((i-1)*2*pi/len)**2*(-Fcoef(2*i-2)*cos((i-1)*2*pi*zz/len)&
                  - Fcoef(2*i-1)*sin((i-1)*2*pi*zz/len))
          ezppp = ezppp+((i-1)*2*pi/len)**3*(Fcoef(2*i-2)*sin((i-1)*2*pi*zz/len)&
                  - Fcoef(2*i-1)*cos((i-1)*2*pi*zz/len))
        enddo
        ez1=ez1*escale
        ezp1=ezp1*escale
        ezpp1=ezpp1*escale
        ezppp = ezppp*escale
!-----------------------------------------------------------------
        tt = pos(4)
        tmpcos = cos(ww*tt+theta0)
        tmpsin = sin(ww*tt+theta0)
        f1 = -(ezpp1+ez1/xl/xl)/4
        f1p = -(ezppp+ezp1/xl/xl)/4
        r2 = tmp(1)**2+tmp(2)**2
        tmpex = -tmp(1)*(ezp1/2+f1p*r2/4)*tmpcos
        tmpey = -tmp(2)*(ezp1/2+f1p*r2/4)*tmpcos
        tmpez = (ez1+f1*r2)*tmpcos
        tmpbx = tmp(2)/(xl*clite)*(ez1/2+f1*r2/4)*tmpsin
        tmpby = -tmp(1)/(xl*clite)*(ez1/2+f1*r2/4)*tmpsin
        tmpbz = 0.0

        !for E field
        tmp(1) = tmpex
        tmp(2) = tmpey*cos(anglex)-tmpez*sin(anglex)
        tmp(3) = tmpey*sin(anglex)+tmpez*cos(anglex)
        temp(1) = tmp(1)*cos(angley)-tmp(3)*sin(angley)
        temp(2) = tmp(2)
        temp(3) = tmp(1)*sin(angley)+tmp(3)*cos(angley)
        extfld(1) = extfld(1)+temp(1)*cos(anglez) - temp(2)*sin(anglez)
        extfld(2) = extfld(2)+temp(1)*sin(anglez) + temp(2)*cos(anglez)
        extfld(3) = extfld(3)+temp(3)

        !for B field
        tmp(1) = tmpbx
        tmp(2) = tmpby*cos(anglex)-tmpbz*sin(anglex)
        tmp(3) = tmpby*sin(anglex)+tmpbz*cos(anglex)
        temp(1) = tmp(1)*cos(angley)-tmp(3)*sin(angley)
        temp(2) = tmp(2)
        temp(3) = tmp(1)*sin(angley)+tmp(3)*cos(angley)
        extfld(4) = extfld(4)+temp(1)*cos(anglez) - temp(2)*sin(anglez)
        extfld(5) = extfld(5)+temp(1)*sin(anglez) + temp(2)*cos(anglez)
        extfld(6) = extfld(6)+temp(3)

        end subroutine getflderrold_DTL

        !get external field with different field and Quad displacement and rotation errors.
        !Here, there 2 displacement errors and rotation errors for Quad
        subroutine  getflderr_DTL(pos,extfld,this,dx,dy,anglex,angley,&
                                  anglez)
        implicit none
        include 'mpif.h'
        double precision, dimension(4), intent(in) :: pos
        type (DTL), intent(in) :: this
        double precision, intent(in) :: dx,dy,anglex,angley,anglez
        double precision, dimension(6), intent(out) :: extfld
        double precision :: zedge,escale,theta0,ww,len,tt,zz,xl
        double precision :: ez1,ezp1,ezpp1,ezppp,f1,f1p,clite,tmpsin,tmpcos,pi
        double precision :: tmpex,tmpey,tmpez,tmpbx,tmpby,tmpbz
        double precision :: bgrad,bgrad1,bgrad2,len1,len2,r2,zmid
        double precision :: dx1,dy1,anglex1,angley1,anglez1
        double precision :: dx0,dy0,anglex0,angley0,anglez0,zz0
        double precision, dimension(3) :: temp,tmp
        integer :: i

        clite = 299792458.e0
        pi = 2*asin(1.0)

        zedge = this%Param(1)
        escale = this%Param(2)
        len = this%Length
        zmid = zedge + len/2
        ww = this%Param(3)/Scfreq
        xl = Scxl/ww  !real frequency has to be used here
        theta0 = this%Param(4)*asin(1.0)/90
        !move into the tank local coordinate.
        zz0=pos(3)-zedge

        len1 = this%Param(7)
        bgrad1 = this%Param(8)
        len2 = len - this%Param(9)
        bgrad2 = this%Param(10)

        if(zz0.lt.len1) then
          bgrad = bgrad1
        else if(zz0.lt.len2) then
          bgrad = 0.0
        else
          bgrad = bgrad2
        endif

        if(zz0.lt.len1) then
          bgrad = bgrad1
          temp(1) = pos(1) - dx
          temp(2) = pos(2) - dy
          tmp(1) = temp(1)*cos(anglez) + temp(2)*sin(anglez)
          tmp(2) = -temp(1)*sin(anglez) + temp(2)*cos(anglez)
          tmp(3) = pos(3) - zedge
          temp(1) = tmp(1)*cos(angley)+tmp(3)*sin(angley)
          temp(2) = tmp(2)
          temp(3) = -tmp(1)*sin(angley)+tmp(3)*cos(angley)
          tmp(1) = temp(1)
          tmp(2) = temp(2)*cos(anglex)+temp(3)*sin(anglex)
          tmp(3) = -temp(2)*sin(anglex)+temp(3)*cos(anglex)
          zz = tmp(3)

          extfld(1) = 0.0
          extfld(2) = 0.0
          extfld(3) = 0.0
          extfld(4) = bgrad*tmp(2)
          extfld(5) = bgrad*tmp(1)
          extfld(6) = 0.0

          tmp(1) = extfld(4)
          tmp(2) = extfld(5)*cos(anglex)-extfld(6)*sin(anglex)
          tmp(3) = extfld(5)*sin(anglex)+extfld(6)*cos(anglex)
          temp(1) = tmp(1)*cos(angley)-tmp(3)*sin(angley)
          temp(2) = tmp(2)
          temp(3) = tmp(1)*sin(angley)+tmp(3)*cos(angley)
          extfld(4) = temp(1)*cos(anglez) - temp(2)*sin(anglez)
          extfld(5) = temp(1)*sin(anglez) + temp(2)*cos(anglez)
          extfld(6) = temp(3)
        else if(zz0.lt.len2) then
          bgrad = 0.0
          extfld = 0.0
        else
          dx0 = this%Param(16)
          dy0 = this%Param(17)
          anglex0 = this%Param(18)
          angley0 = this%Param(19)
          anglez0 = this%Param(20)
          temp(1) = pos(1) - dx0
          temp(2) = pos(2) - dy0
          tmp(1) = temp(1)*cos(anglez0) + temp(2)*sin(anglez0)
          tmp(2) = -temp(1)*sin(anglez0) + temp(2)*cos(anglez0)
          tmp(3) = pos(3) - zedge
          temp(1) = tmp(1)*cos(angley0)+tmp(3)*sin(angley0)
          temp(2) = tmp(2)
          temp(3) = -tmp(1)*sin(angley0)+tmp(3)*cos(angley0)
          tmp(1) = temp(1)
          tmp(2) = temp(2)*cos(anglex0)+temp(3)*sin(anglex0)
          tmp(3) = -temp(2)*sin(anglex0)+temp(3)*cos(anglex0)
          zz = tmp(3)

          extfld(1) = 0.0
          extfld(2) = 0.0
          extfld(3) = 0.0
          extfld(4) = bgrad*tmp(2)
          extfld(5) = bgrad*tmp(1)
          extfld(6) = 0.0

          tmp(1) = extfld(4)
          tmp(2) = extfld(5)*cos(anglex0)-extfld(6)*sin(anglex0)
          tmp(3) = extfld(5)*sin(anglex0)+extfld(6)*cos(anglex0)
          temp(1) = tmp(1)*cos(angley0)-tmp(3)*sin(angley0)
          temp(2) = tmp(2)
          temp(3) = tmp(1)*sin(angley0)+tmp(3)*cos(angley0)
          extfld(4) = temp(1)*cos(anglez0) - temp(2)*sin(anglez0)
          extfld(5) = temp(1)*sin(anglez0) + temp(2)*cos(anglez0)
          extfld(6) = temp(3)
        endif
! get external RF field on axis from analytical function
!-----------------------------------------------------------------
        dx1 = this%Param(21)
        dy1 = this%Param(22)
        anglex1 = this%Param(23)
        angley1 = this%Param(24)
        anglez1 = this%Param(25)

        zz=pos(3)-zmid
        temp(1) = pos(1) - dx
        temp(2) = pos(2) - dy
        tmp(1) = temp(1)*cos(anglez) + temp(2)*sin(anglez)
        tmp(2) = -temp(1)*sin(anglez) + temp(2)*cos(anglez)
        tmp(3) = zz
        temp(1) = tmp(1)*cos(angley)+tmp(3)*sin(angley)
        temp(2) = tmp(2)
        temp(3) = -tmp(1)*sin(angley)+tmp(3)*cos(angley)
        tmp(1) = temp(1)
        tmp(2) = temp(2)*cos(anglex)+temp(3)*sin(anglex)
        tmp(3) = -temp(2)*sin(anglex)+temp(3)*cos(anglex)
        zz = tmp(3)

        ez1 = Fcoef(1)/2
        ezp1 = 0.0
        ezpp1 = 0.0
        ezppp = 0.0
        do i = 2, (Ndata-1)/2 + 1
          ez1 = ez1 + Fcoef(2*i-2)*cos((i-1)*2*pi*zz/len) + &
                Fcoef(2*i-1)*sin((i-1)*2*pi*zz/len)
          ezp1 = ezp1 + (i-1)*2*pi/len*(-Fcoef(2*i-2)*sin((i-1)*2*pi*zz/len)+&
                Fcoef(2*i-1)*cos((i-1)*2*pi*zz/len))
          ezpp1 = ezpp1+((i-1)*2*pi/len)**2*(-Fcoef(2*i-2)*cos((i-1)*2*pi*zz/len)&
                  - Fcoef(2*i-1)*sin((i-1)*2*pi*zz/len))
          ezppp = ezppp+((i-1)*2*pi/len)**3*(Fcoef(2*i-2)*sin((i-1)*2*pi*zz/len)&
                  - Fcoef(2*i-1)*cos((i-1)*2*pi*zz/len))
        enddo
        ez1=ez1*escale
        ezp1=ezp1*escale
        ezpp1=ezpp1*escale
        ezppp = ezppp*escale
!-----------------------------------------------------------------
        tt = pos(4)
        tmpcos = cos(ww*tt+theta0)
        tmpsin = sin(ww*tt+theta0)
        f1 = -(ezpp1+ez1/xl/xl)/4
        f1p = -(ezppp+ezp1/xl/xl)/4
        r2 = tmp(1)**2+tmp(2)**2
        tmpex = -tmp(1)*(ezp1/2+f1p*r2/4)*tmpcos
        tmpey = -tmp(2)*(ezp1/2+f1p*r2/4)*tmpcos
        tmpez = (ez1+f1*r2)*tmpcos
        tmpbx = tmp(2)/(xl*clite)*(ez1/2+f1*r2/4)*tmpsin
        tmpby = -tmp(1)/(xl*clite)*(ez1/2+f1*r2/4)*tmpsin
        tmpbz = 0.0

        !for E field
        tmp(1) = tmpex
        tmp(2) = tmpey*cos(anglex)-tmpez*sin(anglex)
        tmp(3) = tmpey*sin(anglex)+tmpez*cos(anglex)
        temp(1) = tmp(1)*cos(angley)-tmp(3)*sin(angley)
        temp(2) = tmp(2)
        temp(3) = tmp(1)*sin(angley)+tmp(3)*cos(angley)
        extfld(1) = extfld(1)+temp(1)*cos(anglez) - temp(2)*sin(anglez)
        extfld(2) = extfld(2)+temp(1)*sin(anglez) + temp(2)*cos(anglez)
        extfld(3) = extfld(3)+temp(3)

        !for B field
        tmp(1) = tmpbx
        tmp(2) = tmpby*cos(anglex)-tmpbz*sin(anglex)
        tmp(3) = tmpby*sin(anglex)+tmpbz*cos(anglex)
        temp(1) = tmp(1)*cos(angley)-tmp(3)*sin(angley)
        temp(2) = tmp(2)
        temp(3) = tmp(1)*sin(angley)+tmp(3)*cos(angley)
        extfld(4) = extfld(4)+temp(1)*cos(anglez) - temp(2)*sin(anglez)
        extfld(5) = extfld(5)+temp(1)*sin(anglez) + temp(2)*cos(anglez)
        extfld(6) = extfld(6)+temp(3)

        end subroutine getflderr_DTL

        !get external field without displacement and rotation errors.
        subroutine  getfld_DTL(pos,extfld,this)
        implicit none
        include 'mpif.h'
        double precision, dimension(4), intent(in) :: pos
        type (DTL), intent(in) :: this
        double precision, dimension(6), intent(out) :: extfld
        double precision :: zedge,escale,theta0,ww,len,tt,zz,xl
        double precision :: ez1,ezp1,ezpp1,ezppp,f1,f1p,clite,tmpsin,tmpcos,pi
        double precision :: tmpex,tmpey,tmpez,tmpbx,tmpby,tmpbz
        double precision:: bgrad,bgrad1,bgrad2,len1,len2,r2,zmid
        integer :: i

        clite = 299792458.e0
        pi = 2*asin(1.0)

        zedge = this%Param(1)
        escale = this%Param(2)
        len = this%Length
        zmid = zedge + len/2
        ww = this%Param(3)/Scfreq
        xl = Scxl/ww  !real frequency has to be used here
        theta0 = this%Param(4)*asin(1.0)/90
        !move into the tank local coordinate.
        zz=pos(3)-zedge
        len1 = this%Param(7)
        bgrad1 = this%Param(8)
        len2 = len - this%Param(9)
        bgrad2 = this%Param(10)
        if(zz.lt.len1) then
          bgrad = bgrad1
        else if(zz.lt.len2) then
          bgrad = 0.0
        else
          bgrad = bgrad2
        endif
        extfld(1) = 0.0
        extfld(2) = 0.0
        extfld(3) = 0.0
        extfld(4) = bgrad*pos(2)
        extfld(5) = bgrad*pos(1)
        extfld(6) = 0.0
        !This search is based on the assumption that the data is uniformly
        !distributed along z.
! get external RF field on axis from analytical function
!-----------------------------------------------------------------
        zz=pos(3)-zmid
        ez1 = Fcoef(1)/2
        ezp1 = 0.0
        ezpp1 = 0.0
        ezppp = 0.0
        do i = 2, (Ndata-1)/2 + 1
          ez1 = ez1 + Fcoef(2*i-2)*cos((i-1)*2*pi*zz/len) + &
                Fcoef(2*i-1)*sin((i-1)*2*pi*zz/len)
          ezp1 = ezp1 + (i-1)*2*pi/len*(-Fcoef(2*i-2)*sin((i-1)*2*pi*zz/len)+&
                Fcoef(2*i-1)*cos((i-1)*2*pi*zz/len))
          ezpp1 = ezpp1+((i-1)*2*pi/len)**2*(-Fcoef(2*i-2)*cos((i-1)*2*pi*zz/len)&
                  - Fcoef(2*i-1)*sin((i-1)*2*pi*zz/len))
          ezppp = ezppp+((i-1)*2*pi/len)**3*(Fcoef(2*i-2)*sin((i-1)*2*pi*zz/len)&
                  - Fcoef(2*i-1)*cos((i-1)*2*pi*zz/len))
        enddo
        ez1=ez1*escale
        ezp1=ezp1*escale
        ezpp1=ezpp1*escale
        ezppp = ezppp*escale
!-----------------------------------------------------------------
        tt = pos(4)
        tmpcos = cos(ww*tt+theta0)
        tmpsin = sin(ww*tt+theta0)
        f1 = -(ezpp1+ez1/xl/xl)/4
        f1p = -(ezppp+ezp1/xl/xl)/4
!        f1 = 0.0
!        f1p = 0.0
        r2 = pos(1)**2+pos(2)**2
        tmpex = -pos(1)*(ezp1/2+f1p*r2/4)*tmpcos
        tmpey = -pos(2)*(ezp1/2+f1p*r2/4)*tmpcos
        tmpez = (ez1+f1*r2)*tmpcos
        tmpbx = pos(2)/(xl*clite)*(ez1/2+f1*r2/4)*tmpsin
        tmpby = -pos(1)/(xl*clite)*(ez1/2+f1*r2/4)*tmpsin
        tmpbz = 0.0
        extfld(1) = extfld(1) + tmpex
        extfld(2) = extfld(2) + tmpey
        extfld(3) = extfld(3) + tmpez
        extfld(4) = extfld(4) + tmpbx
        extfld(5) = extfld(5) + tmpby
        extfld(6) = extfld(6) + tmpbz

        end subroutine getfld_DTL

        !get external field without displacement and rotation errors.
        subroutine  getfldpts_DTL(pos,extfld,this,nplc)
        implicit none
        include 'mpif.h'
        integer,intent(in) ::nplc
        double precision, dimension(4,nplc), intent(in) :: pos
        type (DTL), intent(in) :: this
        double precision, dimension(6,nplc), intent(out) :: extfld
        double precision :: zedge,escale,theta0,ww,len,tt,zz,xl
        double precision :: ez1,ezp1,ezpp1,ezppp,f1,f1p,clite,tmpsin,tmpcos,pi
        double precision :: tmpex,tmpey,tmpez,tmpbx,tmpby,tmpbz
        double precision:: bgrad,bgrad1,bgrad2,len1,len2,r2,zmid
        integer :: i,ipt

        clite = 299792458.e0
        pi = 2*asin(1.0)

        zedge = this%Param(1)
        escale = this%Param(2)
        len = this%Length
        zmid = zedge + len/2
        ww = this%Param(3)/Scfreq
        xl = Scxl/ww  !real frequency has to be used here
        theta0 = this%Param(4)*asin(1.0)/90
        len1 = this%Param(7)
        bgrad1 = this%Param(8)
        len2 = len - this%Param(9)
        bgrad2 = this%Param(10)

        !move into the tank local coordinate.
        zz=pos(3,1)-zedge
        if(zz.lt.len1) then
          bgrad = bgrad1
        else if(zz.lt.len2) then
          bgrad = 0.0
        else
          bgrad = bgrad2
        endif

        do ipt =1,nplc

        extfld(1,ipt) = 0.0
        extfld(2,ipt) = 0.0
        extfld(3,ipt) = 0.0
        extfld(4,ipt) = bgrad*pos(2,ipt)
        extfld(5,ipt) = bgrad*pos(1,ipt)
        extfld(6,ipt) = 0.0

        end do

        !This search is based on the assumption that the data is uniformly
        !distributed along z.
! get external RF field on axis from analytical function
!-----------------------------------------------------------------
        zz=pos(3,1)-zmid
        ez1 = Fcoef(1)/2
        ezp1 = 0.0
        ezpp1 = 0.0
        ezppp = 0.0
        do i = 2, (Ndata-1)/2 + 1
          ez1 = ez1 + Fcoef(2*i-2)*cos((i-1)*2*pi*zz/len) + &
                Fcoef(2*i-1)*sin((i-1)*2*pi*zz/len)
          ezp1 = ezp1 + (i-1)*2*pi/len*(-Fcoef(2*i-2)*sin((i-1)*2*pi*zz/len)+&
                Fcoef(2*i-1)*cos((i-1)*2*pi*zz/len))
          ezpp1 = ezpp1+((i-1)*2*pi/len)**2*(-Fcoef(2*i-2)*cos((i-1)*2*pi*zz/len)&
                  - Fcoef(2*i-1)*sin((i-1)*2*pi*zz/len))
          ezppp = ezppp+((i-1)*2*pi/len)**3*(Fcoef(2*i-2)*sin((i-1)*2*pi*zz/len)&
                  - Fcoef(2*i-1)*cos((i-1)*2*pi*zz/len))
        enddo
        ez1=ez1*escale
        ezp1=ezp1*escale
        ezpp1=ezpp1*escale
        ezppp = ezppp*escale
!-----------------------------------------------------------------
        f1 = -(ezpp1+ez1/xl/xl)/4
        f1p = -(ezppp+ezp1/xl/xl)/4
!        f1 = 0.0
!        f1p = 0.0

        do ipt=1,nplc

        tt = pos(4,ipt)
        tmpcos = cos(ww*tt+theta0)
        tmpsin = sin(ww*tt+theta0)
        r2 = pos(1,ipt)**2+pos(2,ipt)**2
        tmpex = -pos(1,ipt)*(ezp1/2+f1p*r2/4)*tmpcos
        tmpey = -pos(2,ipt)*(ezp1/2+f1p*r2/4)*tmpcos
        tmpez = (ez1+f1*r2)*tmpcos
        tmpbx = pos(2,ipt)/(xl*clite)*(ez1/2+f1*r2/4)*tmpsin
        tmpby = -pos(1,ipt)/(xl*clite)*(ez1/2+f1*r2/4)*tmpsin
        tmpbz = 0.0
        extfld(1,ipt) = extfld(1,ipt) + tmpex
        extfld(2,ipt) = extfld(2,ipt) + tmpey
        extfld(3,ipt) = extfld(3,ipt) + tmpez
        extfld(4,ipt) = extfld(4,ipt) + tmpbx
        extfld(5,ipt) = extfld(5,ipt) + tmpby
        extfld(6,ipt) = extfld(6,ipt) + tmpbz

        enddo

        end subroutine getfldpts_DTL

        !get external field with different field and Quad displacement and rotation errors.
        !Here, there 2 displacement errors and rotation errors for Quad
        subroutine  getflderrpts_DTL(pos,extfld,this,dx,dy,anglex,angley,&
                                  anglez,nplc)
        implicit none
        include 'mpif.h'
        integer,intent(in)::nplc
        double precision, dimension(4,nplc), intent(in) :: pos
        type (DTL), intent(in) :: this
        double precision, intent(in) :: dx,dy,anglex,angley,anglez
        double precision, dimension(6,nplc), intent(out) :: extfld
        double precision :: zedge,escale,theta0,ww,len,tt,zz,xl
        double precision :: ez1,ezp1,ezpp1,ezppp,f1,f1p,clite,tmpsin,tmpcos,pi
        double precision :: tmpex,tmpey,tmpez,tmpbx,tmpby,tmpbz
        double precision :: bgrad,bgrad1,bgrad2,len1,len2,r2,zmid
        double precision :: dx1,dy1,anglex1,angley1,anglez1
        double precision :: dx0,dy0,anglex0,angley0,anglez0,zz0
        double precision, dimension(3) :: temp,tmp
        integer :: i,ipt

        clite = 299792458.e0
        pi = 2*asin(1.0)

        zedge = this%Param(1)
        escale = this%Param(2)
        len = this%Length
        zmid = zedge + len/2
        ww = this%Param(3)/Scfreq
        xl = Scxl/ww  !real frequency has to be used here
        theta0 = this%Param(4)*asin(1.0)/90
        len1 = this%Param(7)
        bgrad1 = this%Param(8)
        len2 = len - this%Param(9)
        bgrad2 = this%Param(10)

        !move into the tank local coordinate.
        zz0=pos(3,1)-zedge

        if(zz0.lt.len1) then
          bgrad = bgrad1
        else if(zz0.lt.len2) then
          bgrad = 0.0
        else
          bgrad = bgrad2
        endif

        do ipt=1,nplc

        if(zz0.lt.len1) then
          bgrad = bgrad1
          temp(1) = pos(1,ipt) - dx
          temp(2) = pos(2,ipt) - dy
          tmp(1) = temp(1)*cos(anglez) + temp(2)*sin(anglez)
          tmp(2) = -temp(1)*sin(anglez) + temp(2)*cos(anglez)
          tmp(3) = pos(3,ipt) - zedge
          temp(1) = tmp(1)*cos(angley)+tmp(3)*sin(angley)
          temp(2) = tmp(2)
          temp(3) = -tmp(1)*sin(angley)+tmp(3)*cos(angley)
          tmp(1) = temp(1)
          tmp(2) = temp(2)*cos(anglex)+temp(3)*sin(anglex)
          tmp(3) = -temp(2)*sin(anglex)+temp(3)*cos(anglex)
          zz = tmp(3)

          extfld(1,ipt) = 0.0
          extfld(2,ipt) = 0.0
          extfld(3,ipt) = 0.0
          extfld(4,ipt) = bgrad*tmp(2)
          extfld(5,ipt) = bgrad*tmp(1)
          extfld(6,ipt) = 0.0

          tmp(1) = extfld(4,ipt)
          tmp(2) = extfld(5,ipt)*cos(anglex)-extfld(6,ipt)*sin(anglex)
          tmp(3) = extfld(5,ipt)*sin(anglex)+extfld(6,ipt)*cos(anglex)
          temp(1) = tmp(1)*cos(angley)-tmp(3)*sin(angley)
          temp(2) = tmp(2)
          temp(3) = tmp(1)*sin(angley)+tmp(3)*cos(angley)
          extfld(4,ipt) = temp(1)*cos(anglez) - temp(2)*sin(anglez)
          extfld(5,ipt) = temp(1)*sin(anglez) + temp(2)*cos(anglez)
          extfld(6,ipt) = temp(3)
        else if(zz0.lt.len2) then
          bgrad = 0.0
          extfld(:,ipt) = 0.0
        else
          dx0 = this%Param(16)
          dy0 = this%Param(17)
          anglex0 = this%Param(18)
          angley0 = this%Param(19)
          anglez0 = this%Param(20)
          temp(1) = pos(1,ipt) - dx0
          temp(2) = pos(2,ipt) - dy0
          tmp(1) = temp(1)*cos(anglez0) + temp(2)*sin(anglez0)
          tmp(2) = -temp(1)*sin(anglez0) + temp(2)*cos(anglez0)
          tmp(3) = pos(3,ipt) - zedge
          temp(1) = tmp(1)*cos(angley0)+tmp(3)*sin(angley0)
          temp(2) = tmp(2)
          temp(3) = -tmp(1)*sin(angley0)+tmp(3)*cos(angley0)
          tmp(1) = temp(1)
          tmp(2) = temp(2)*cos(anglex0)+temp(3)*sin(anglex0)
          tmp(3) = -temp(2)*sin(anglex0)+temp(3)*cos(anglex0)
          zz = tmp(3)

          extfld(1,ipt) = 0.0
          extfld(2,ipt) = 0.0
          extfld(3,ipt) = 0.0
          extfld(4,ipt) = bgrad*tmp(2)
          extfld(5,ipt) = bgrad*tmp(1)
          extfld(6,ipt) = 0.0

          tmp(1) = extfld(4,ipt)
          tmp(2) = extfld(5,ipt)*cos(anglex0)-extfld(6,ipt)*sin(anglex0)
          tmp(3) = extfld(5,ipt)*sin(anglex0)+extfld(6,ipt)*cos(anglex0)
          temp(1) = tmp(1)*cos(angley0)-tmp(3)*sin(angley0)
          temp(2) = tmp(2)
          temp(3) = tmp(1)*sin(angley0)+tmp(3)*cos(angley0)
          extfld(4,ipt) = temp(1)*cos(anglez0) - temp(2)*sin(anglez0)
          extfld(5,ipt) = temp(1)*sin(anglez0) + temp(2)*cos(anglez0)
          extfld(6,ipt) = temp(3)
        endif
! get external RF field on axis from analytical function
!-----------------------------------------------------------------
        dx1 = this%Param(21)
        dy1 = this%Param(22)
        anglex1 = this%Param(23)
        angley1 = this%Param(24)
        anglez1 = this%Param(25)

        zz=pos(3,ipt)-zmid
        temp(1) = pos(1,ipt) - dx
        temp(2) = pos(2,ipt) - dy
        tmp(1) = temp(1)*cos(anglez) + temp(2)*sin(anglez)
        tmp(2) = -temp(1)*sin(anglez) + temp(2)*cos(anglez)
        tmp(3) = zz
        temp(1) = tmp(1)*cos(angley)+tmp(3)*sin(angley)
        temp(2) = tmp(2)
        temp(3) = -tmp(1)*sin(angley)+tmp(3)*cos(angley)
        tmp(1) = temp(1)
        tmp(2) = temp(2)*cos(anglex)+temp(3)*sin(anglex)
        tmp(3) = -temp(2)*sin(anglex)+temp(3)*cos(anglex)
        zz = tmp(3)

        ez1 = Fcoef(1)/2
        ezp1 = 0.0
        ezpp1 = 0.0
        ezppp = 0.0
        do i = 2, (Ndata-1)/2 + 1
          ez1 = ez1 + Fcoef(2*i-2)*cos((i-1)*2*pi*zz/len) + &
                Fcoef(2*i-1)*sin((i-1)*2*pi*zz/len)
          ezp1 = ezp1 + (i-1)*2*pi/len*(-Fcoef(2*i-2)*sin((i-1)*2*pi*zz/len)+&
                Fcoef(2*i-1)*cos((i-1)*2*pi*zz/len))
          ezpp1 = ezpp1+((i-1)*2*pi/len)**2*(-Fcoef(2*i-2)*cos((i-1)*2*pi*zz/len)&
                  - Fcoef(2*i-1)*sin((i-1)*2*pi*zz/len))
          ezppp = ezppp+((i-1)*2*pi/len)**3*(Fcoef(2*i-2)*sin((i-1)*2*pi*zz/len)&
                  - Fcoef(2*i-1)*cos((i-1)*2*pi*zz/len))
        enddo
        ez1=ez1*escale
        ezp1=ezp1*escale
        ezpp1=ezpp1*escale
        ezppp = ezppp*escale
!-----------------------------------------------------------------
        tt = pos(4,ipt)
        tmpcos = cos(ww*tt+theta0)
        tmpsin = sin(ww*tt+theta0)
        f1 = -(ezpp1+ez1/xl/xl)/4
        f1p = -(ezppp+ezp1/xl/xl)/4
        r2 = tmp(1)**2+tmp(2)**2
        tmpex = -tmp(1)*(ezp1/2+f1p*r2/4)*tmpcos
        tmpey = -tmp(2)*(ezp1/2+f1p*r2/4)*tmpcos
        tmpez = (ez1+f1*r2)*tmpcos
        tmpbx = tmp(2)/(xl*clite)*(ez1/2+f1*r2/4)*tmpsin
        tmpby = -tmp(1)/(xl*clite)*(ez1/2+f1*r2/4)*tmpsin
        tmpbz = 0.0

        !for E field
        tmp(1) = tmpex
        tmp(2) = tmpey*cos(anglex)-tmpez*sin(anglex)
        tmp(3) = tmpey*sin(anglex)+tmpez*cos(anglex)
        temp(1) = tmp(1)*cos(angley)-tmp(3)*sin(angley)
        temp(2) = tmp(2)
        temp(3) = tmp(1)*sin(angley)+tmp(3)*cos(angley)
        extfld(1,ipt) = extfld(1,ipt)+temp(1)*cos(anglez) - temp(2)*sin(anglez)
        extfld(2,ipt) = extfld(2,ipt)+temp(1)*sin(anglez) + temp(2)*cos(anglez)
        extfld(3,ipt) = extfld(3,ipt)+temp(3)

        !for B field
        tmp(1) = tmpbx
        tmp(2) = tmpby*cos(anglex)-tmpbz*sin(anglex)
        tmp(3) = tmpby*sin(anglex)+tmpbz*cos(anglex)
        temp(1) = tmp(1)*cos(angley)-tmp(3)*sin(angley)
        temp(2) = tmp(2)
        temp(3) = tmp(1)*sin(angley)+tmp(3)*cos(angley)
        extfld(4,ipt) = extfld(4,ipt)+temp(1)*cos(anglez) - temp(2)*sin(anglez)
        extfld(5,ipt) = extfld(5,ipt)+temp(1)*sin(anglez) + temp(2)*cos(anglez)
        extfld(6,ipt) = extfld(6,ipt)+temp(3)

        enddo

        end subroutine getflderrpts_DTL
      end module DTLclass

max_stars_repo_name: hhslepicka/IMPACT-Z

max_issues_repo_name: hhslepicka/IMPACT-Z

max_forks_repo_name: hhslepicka/IMPACT-Z


program tarefab1

    ! Número de passos N
    parameter (N = 1000)
    ! vetor posição
    dimension iposicao(-N:N)

    ! definindo arquivo de saida
    open(10, file='saida-b1-10407962')

    ! recebe a probabilidade p de ir para a direita
    write(*,*) "valor p:"
    read(*,*) p

    ! recebe o número de andarilhos M
    print *, "número de andarilhos M:"
    read(*,*) M

    ! defini inicialmente como zero as variáveis
    iposicao = 0
    iandarilho = 0
    soma = 0
    soma2 = 0

    ! loop dos andarilhos
    do j = 1, M
        ! zera a posição do andarilho
        ix = 0

        ! loop dos passos
        do i = 1, N

            ! condicional para direção do passo
            if (rand().lt.p) then
                ! <= p, vai para direita
                ix = ix + 1
            else
                ! > p, vai para esquerda
                ix = ix - 1
            end if
        end do

        ! soma 1 andarilhos na posição que parou
        iposicao(ix) = iposicao(ix) + 1

        ! soma para <x>
        soma = soma + ix
        ! soma2 para <x**2>
        soma2 = soma2 + ix**2
    end do

    ! loop para salvar as posições no arquivo de saida
    do i = -N, N

        ! condicional para o caso de não tiver andarilho na posição
        ! não salva no arquivo
        if ( iposicao(i).ne.0 ) then
            write(10,'(I0," ",I0)') i, iposicao(i)
        end if
    end do

    ! escreve no terminal o número de passos, <x> e <x**2>
    write(*,'("N: ",I0,"  <x>: ",1F0.3," <x**2>: ",1F0.3)') N, abs(soma)/M, abs(soma2)/M

    ! fecha o arquivo de saida
    close(10)

 end program tarefab1


PROGRAM non_block

  USE mpi
  IMPLICIT NONE

  INTEGER, PARAMETER :: tag = 100
  INTEGER, PARAMETER :: n_elements = 1000

  INTEGER :: rank
  INTEGER, DIMENSION(n_elements, n_elements) :: send_val, recv_val
  INTEGER :: nproc
  INTEGER :: left, right
  INTEGER, DIMENSION(2) :: requests
  INTEGER, DIMENSION(MPI_STATUS_SIZE, 2) :: statuses
  INTEGER :: ierr
  REAL :: time1, time2, dtime

  CALL MPI_Init(ierr)

  CALL MPI_Comm_size(MPI_COMM_WORLD, nproc, ierr)
  CALL MPI_Comm_rank(MPI_COMM_WORLD, rank, ierr)

  !Set up periodic domain
  left = rank - 1
  IF (left < 0) left = nproc - 1
  right = rank + 1
  IF (right > nproc - 1) right = 0

  send_val = rank
  CALL MPI_Isend(send_val(1,:), n_elements, MPI_INTEGER, right, tag, &
      MPI_COMM_WORLD, requests(1), ierr)
  CALL MPI_Irecv(recv_val(1,:), n_elements, MPI_INTEGER, left, tag, &
      MPI_COMM_WORLD, requests(2), ierr)

  CALL MPI_Waitall(2, requests, statuses, ierr)

  PRINT *,"Rank ", rank, " got minimal value from rank ", left, " of ", &
      MINVAL(recv_val(1,:))
  CALL MPI_Barrier(MPI_COMM_WORLD, ierr)
  IF (rank == 0) PRINT *, "It took ", dtime, " seconds to complete"

  CALL MPI_Finalize(ierr)

END PROGRAM non_block


max_stars_repo_name: umarcor/svunit

max_issues_repo_name: umarcor/svunit

max_forks_repo_name: umarcor/svunit


text:
! The MIT License (MIT)
!
! Copyright (c) 2016 odlomax
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.

! string manipulation module
module m_string

   use m_kind_parameters

   implicit none

   contains

   ! split space delineated string up into reals
   pure function string_to_real(string) result (value)

      ! argument declarations
      character(kind=chr_kind,len=string_length),intent(in) :: string     ! character string

      ! result declaration
      real(kind=rel_kind),allocatable :: value(:)                          ! array of reals

      ! variable declarations
      integer(kind=int_kind) :: j                                          ! counter
      integer(kind=int_kind) :: i_sub_string                               ! position of sub_string in string
      integer(kind=int_kind) :: n_value                                    ! number of entries in string
      character(kind=chr_kind,len=string_length) :: trimmed_string         ! trimmed version of string variable

      ! count number of sub_string delineated entries
      trimmed_string=adjustl(string)
      n_value=0
      do

         if (len(trim(adjustl(trimmed_string)))==0) exit

         i_sub_string=index(trimmed_string," ")
         n_value=n_value+1
         trimmed_string=adjustl(trimmed_string(i_sub_string+1:))

      end do

      ! read string and assign to value array
      allocate(value(n_value))
      trimmed_string=adjustl(string)
      do j=1,n_value

         i_sub_string=index(trimmed_string," ")
         read(trimmed_string(:i_sub_string-1),*) value(j)
         trimmed_string=adjustl(trimmed_string(i_sub_string+1:))

      end do

      return

   end function

end module

max_stars_repo_name: odlomax/spamcart-dev

max_issues_repo_name: odlomax/spamcart-dev

max_forks_repo_name: odlomax/spamcart-dev


text: #redirect wiki:wikispot:Help_with_Editing


text:
c Convert latitude and longitude to kilometers relative
c to center of coordinates by short distance conversion.

	subroutine dist(xlat, xlon, xkm, ykm)

	implicit none

c	Parameters:
	doubleprecision	xlat, xlon	! (input)
	real		xkm, ykm	! (output)

c	Local variables:
	doubleprecision lat1, lat2, lat3
	real	q
	real	xx
	real	yp

	include "geocoord.inc"

c Set up short distance conversion by subr. SETORG
      q=60*xlat-olat
      yp=q+olat
      lat1=datan(rlatc*dtan(RAD*yp/60.0))
      lat2=datan(rlatc*dtan(RAD*OLAT/60.0))
      LAT3=(LAT2+LAT1)/2.
      xx=60*xlon-olon  !  - wegen LON E
      q=q*aa
      xx = xx*bb*dcos(LAT3)
      IF(rotate.ne.0.) then
c** rotate coordinate system anticlockwise
        yp=cost*q+sint*xx
        xx=cost*xx-sint*q
        q=yp
      ENDIF

      xkm=xx
      ykm=q

      return
      end