X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=blobdiff_plain;f=source%2Funres%2Fsrc_CSA_DiL%2Flocal_move.f;fp=source%2Funres%2Fsrc_CSA_DiL%2Flocal_move.f;h=0000000000000000000000000000000000000000;hp=763d3cc821dfd76db68258f1b1a3ee0ba19acb3e;hb=de4bc5453ea46e111d936cb85e1758ed21c08fcd;hpb=b75425747e3e2b448ca5e0ef8367712e1f339124 diff --git a/source/unres/src_CSA_DiL/local_move.f b/source/unres/src_CSA_DiL/local_move.f deleted file mode 100644 index 763d3cc..0000000 --- a/source/unres/src_CSA_DiL/local_move.f +++ /dev/null @@ -1,972 +0,0 @@ -c------------------------------------------------------------- - - subroutine local_move_init(debug) -crc implicit none - -c Includes - implicit real*8 (a-h,o-z) - include 'DIMENSIONS' ! Needed by COMMON.LOCAL - include 'COMMON.GEO' ! For pi, deg2rad - include 'COMMON.LOCAL' ! For vbl - include 'COMMON.LOCMOVE' - -c INPUT arguments - logical debug - - -c Determine wheter to do some debugging output - locmove_output=debug - -c Set the init_called flag to 1 - init_called=1 - -c The following are never changed - min_theta=60.D0*deg2rad ! (0,PI) - max_theta=175.D0*deg2rad ! (0,PI) - dmin2=vbl*vbl*2.*(1.-cos(min_theta)) - dmax2=vbl*vbl*2.*(1.-cos(max_theta)) - flag=1.0D300 - small=1.0D-5 - small2=0.5*small*small - -c Not really necessary... - a_n=0 - b_n=0 - res_n=0 - - return - end - -c------------------------------------------------------------- - - subroutine local_move(n_start, n_end, PHImin, PHImax) -c Perform a local move between residues m and n (inclusive) -c PHImin and PHImax [0,PI] determine the size of the move -c Works on whatever structure is in the variables theta and phi, -c sidechain variables are left untouched -c The final structure is NOT minimized, but both the cartesian -c variables c and the angles are up-to-date at the end (no further -c chainbuild is required) -crc implicit none - -c Includes - implicit real*8 (a-h,o-z) - include 'DIMENSIONS' - include 'COMMON.GEO' - include 'COMMON.CHAIN' - include 'COMMON.VAR' - include 'COMMON.MINIM' - include 'COMMON.SBRIDGE' - include 'COMMON.LOCMOVE' - -c External functions - integer move_res - external move_res - double precision ran_number - external ran_number - -c INPUT arguments - integer n_start, n_end ! First and last residues to move - double precision PHImin, PHImax ! min/max angles [0,PI] - -c Local variables - integer i,j - double precision min,max - integer iretcode - - -c Check if local_move_init was called. This assumes that it -c would not be 1 if not explicitely initialized - if (init_called.ne.1) then - write(6,*)' *** local_move_init not called!!!' - stop - endif - -c Quick check for crazy range - if (n_start.gt.n_end .or. n_start.lt.1 .or. n_end.gt.nres) then - write(6,'(a,i3,a,i3)') - + ' *** Cannot make local move between n_start = ', - + n_start,' and n_end = ',n_end - return - endif - -c Take care of end residues first... - if (n_start.eq.1) then -c Move residue 1 (completely random) - theta(3)=ran_number(min_theta,max_theta) - phi(4)=ran_number(-PI,PI) - i=2 - else - i=n_start - endif - if (n_end.eq.nres) then -c Move residue nres (completely random) - theta(nres)=ran_number(min_theta,max_theta) - phi(nres)=ran_number(-PI,PI) - j=nres-1 - else - j=n_end - endif - -c ...then go through all other residues one by one -c Start from the two extremes and converge - call chainbuild - do while (i.le.j) - min=PHImin - max=PHImax -c$$$c Move the first two residues by less than the others -c$$$ if (i-n_start.lt.3) then -c$$$ if (i-n_start.eq.0) then -c$$$ min=0.4*PHImin -c$$$ max=0.4*PHImax -c$$$ else if (i-n_start.eq.1) then -c$$$ min=0.8*PHImin -c$$$ max=0.8*PHImax -c$$$ else if (i-n_start.eq.2) then -c$$$ min=PHImin -c$$$ max=PHImax -c$$$ endif -c$$$ endif - -c The actual move, on residue i - iretcode=move_res(min,max,i) ! Discard iretcode - i=i+1 - - if (i.le.j) then - min=PHImin - max=PHImax -c$$$c Move the last two residues by less than the others -c$$$ if (n_end-j.lt.3) then -c$$$ if (n_end-j.eq.0) then -c$$$ min=0.4*PHImin -c$$$ max=0.4*PHImax -c$$$ else if (n_end-j.eq.1) then -c$$$ min=0.8*PHImin -c$$$ max=0.8*PHImax -c$$$ else if (n_end-j.eq.2) then -c$$$ min=PHImin -c$$$ max=PHImax -c$$$ endif -c$$$ endif - -c The actual move, on residue j - iretcode=move_res(min,max,j) ! Discard iretcode - j=j-1 - endif - enddo - - call int_from_cart(.false.,.false.) - - return - end - -c------------------------------------------------------------- - - subroutine output_tabs -c Prints out the contents of a_..., b_..., res_... - implicit none - -c Includes - include 'COMMON.GEO' - include 'COMMON.LOCMOVE' - -c Local variables - integer i,j - - - write(6,*)'a_...' - write(6,'(8f7.1)')(a_ang(i)*rad2deg,i=0,a_n-1) - write(6,'(8(2x,3l1,2x))')((a_tab(i,j),i=0,2),j=0,a_n-1) - - write(6,*)'b_...' - write(6,'(4f7.1)')(b_ang(i)*rad2deg,i=0,b_n-1) - write(6,'(4(2x,3l1,2x))')((b_tab(i,j),i=0,2),j=0,b_n-1) - - write(6,*)'res_...' - write(6,'(12f7.1)')(res_ang(i)*rad2deg,i=0,res_n-1) - write(6,'(12(2x,3l1,2x))')((res_tab(0,i,j),i=0,2),j=0,res_n-1) - write(6,'(12(2x,3l1,2x))')((res_tab(1,i,j),i=0,2),j=0,res_n-1) - write(6,'(12(2x,3l1,2x))')((res_tab(2,i,j),i=0,2),j=0,res_n-1) - - return - end - -c------------------------------------------------------------- - - subroutine angles2tab(PHImin,PHImax,n,ang,tab) -c Only uses angles if [0,PI] (but PHImin cannot be 0., -c and PHImax cannot be PI) - implicit none - -c Includes - include 'COMMON.GEO' - -c INPUT arguments - double precision PHImin,PHImax - -c OUTPUT arguments - integer n - double precision ang(0:3) - logical tab(0:2,0:3) - - - if (PHImin .eq. PHImax) then -c Special case with two 010's - n = 2; - ang(0) = -PHImin; - ang(1) = PHImin; - tab(0,0) = .false. - tab(2,0) = .false. - tab(0,1) = .false. - tab(2,1) = .false. - tab(1,0) = .true. - tab(1,1) = .true. - else if (PHImin .eq. PI) then -c Special case with one 010 - n = 1 - ang(0) = PI - tab(0,0) = .false. - tab(2,0) = .false. - tab(1,0) = .true. - else if (PHImax .eq. 0.) then -c Special case with one 010 - n = 1 - ang(0) = 0. - tab(0,0) = .false. - tab(2,0) = .false. - tab(1,0) = .true. - else -c Standard cases - n = 0 - if (PHImin .gt. 0.) then -c Start of range (011) - ang(n) = PHImin - tab(0,n) = .false. - tab(1,n) = .true. - tab(2,n) = .true. -c End of range (110) - ang(n+1) = -PHImin - tab(0,n+1) = .true. - tab(1,n+1) = .true. - tab(2,n+1) = .false. - n = n+2 - endif - if (PHImax .lt. PI) then -c Start of range (011) - ang(n) = -PHImax - tab(0,n) = .false. - tab(1,n) = .true. - tab(2,n) = .true. -c End of range (110) - ang(n+1) = PHImax - tab(0,n+1) = .true. - tab(1,n+1) = .true. - tab(2,n+1) = .false. - n = n+2 - endif - endif - - return - end - -c------------------------------------------------------------- - - subroutine minmax_angles(x,y,z,r,n,ang,tab) -c When solutions do not exist, assume all angles -c are acceptable - i.e., initial geometry must be correct - implicit none - -c Includes - include 'COMMON.GEO' - include 'COMMON.LOCMOVE' - -c Input arguments - double precision x,y,z,r - -c Output arguments - integer n - double precision ang(0:3) - logical tab(0:2,0:3) - -c Local variables - double precision num, denom, phi - double precision Kmin, Kmax - integer i - - - num = x*x + y*y + z*z - denom = x*x + y*y - n = 0 - if (denom .gt. 0.) then - phi = atan2(y,x) - denom = 2.*r*sqrt(denom) - num = num+r*r - Kmin = (num - dmin2)/denom - Kmax = (num - dmax2)/denom - -c Allowed values of K (else all angles are acceptable) -c -1 <= Kmin < 1 -c -1 < Kmax <= 1 - if (Kmin .gt. 1. .or. abs(Kmin-1.) .lt. small2) then - Kmin = -flag - else if (Kmin .lt. -1. .or. abs(Kmin+1.) .lt. small2) then - Kmin = PI - else - Kmin = acos(Kmin) - endif - - if (Kmax .lt. -1. .or. abs(Kmax+1.) .lt. small2) then - Kmax = flag - else if (Kmax .gt. 1. .or. abs(Kmax-1.) .lt. small2) then - Kmax = 0. - else - Kmax = acos(Kmax) - endif - - if (Kmax .lt. Kmin) Kmax = Kmin - - call angles2tab(Kmin, Kmax, n, ang, tab) - -c Add phi and check that angles are within range (-PI,PI] - do i=0,n-1 - ang(i) = ang(i)+phi - if (ang(i) .le. -PI) then - ang(i) = ang(i)+2.*PI - else if (ang(i) .gt. PI) then - ang(i) = ang(i)-2.*PI - endif - enddo - endif - - return - end - -c------------------------------------------------------------- - - subroutine construct_tab -c Take a_... and b_... values and produces the results res_... -c x_ang are assumed to be all different (diff > small) -c x_tab(1,i) must be 1 for all i (i.e., all x_ang are acceptable) - implicit none - -c Includes - include 'COMMON.LOCMOVE' - -c Local variables - integer n_max,i,j,index - logical done - double precision phi - - - n_max = a_n + b_n - if (n_max .eq. 0) then - res_n = 0 - return - endif - - do i=0,n_max-1 - do j=0,1 - res_tab(j,0,i) = .true. - res_tab(j,2,i) = .true. - res_tab(j,1,i) = .false. - enddo - enddo - - index = 0 - phi = -flag - done = .false. - do while (.not.done) - res_ang(index) = flag - -c Check a first... - do i=0,a_n-1 - if ((a_ang(i)-phi).gt.small .and. - + a_ang(i) .lt. res_ang(index)) then -c Found a lower angle - res_ang(index) = a_ang(i) -c Copy the values from a_tab into res_tab(0,,) - res_tab(0,0,index) = a_tab(0,i) - res_tab(0,1,index) = a_tab(1,i) - res_tab(0,2,index) = a_tab(2,i) -c Set default values for res_tab(1,,) - res_tab(1,0,index) = .true. - res_tab(1,1,index) = .false. - res_tab(1,2,index) = .true. - else if (abs(a_ang(i)-res_ang(index)).lt.small) then -c Found an equal angle (can only be equal to a b_ang) - res_tab(0,0,index) = a_tab(0,i) - res_tab(0,1,index) = a_tab(1,i) - res_tab(0,2,index) = a_tab(2,i) - endif - enddo -c ...then check b - do i=0,b_n-1 - if ((b_ang(i)-phi).gt.small .and. - + b_ang(i) .lt. res_ang(index)) then -c Found a lower angle - res_ang(index) = b_ang(i) -c Copy the values from b_tab into res_tab(1,,) - res_tab(1,0,index) = b_tab(0,i) - res_tab(1,1,index) = b_tab(1,i) - res_tab(1,2,index) = b_tab(2,i) -c Set default values for res_tab(0,,) - res_tab(0,0,index) = .true. - res_tab(0,1,index) = .false. - res_tab(0,2,index) = .true. - else if (abs(b_ang(i)-res_ang(index)).lt.small) then -c Found an equal angle (can only be equal to an a_ang) - res_tab(1,0,index) = b_tab(0,i) - res_tab(1,1,index) = b_tab(1,i) - res_tab(1,2,index) = b_tab(2,i) - endif - enddo - - if (res_ang(index) .eq. flag) then - res_n = index - done = .true. - else if (index .eq. n_max-1) then - res_n = n_max - done = .true. - else - phi = res_ang(index) ! Store previous angle - index = index+1 - endif - enddo - -c Fill the gaps -c First a... - index = 0 - if (a_n .gt. 0) then - do while (.not.res_tab(0,1,index)) - index=index+1 - enddo - done = res_tab(0,2,index) - do i=index+1,res_n-1 - if (res_tab(0,1,i)) then - done = res_tab(0,2,i) - else - res_tab(0,0,i) = done - res_tab(0,1,i) = done - res_tab(0,2,i) = done - endif - enddo - done = res_tab(0,0,index) - do i=index-1,0,-1 - if (res_tab(0,1,i)) then - done = res_tab(0,0,i) - else - res_tab(0,0,i) = done - res_tab(0,1,i) = done - res_tab(0,2,i) = done - endif - enddo - else - do i=0,res_n-1 - res_tab(0,0,i) = .true. - res_tab(0,1,i) = .true. - res_tab(0,2,i) = .true. - enddo - endif -c ...then b - index = 0 - if (b_n .gt. 0) then - do while (.not.res_tab(1,1,index)) - index=index+1 - enddo - done = res_tab(1,2,index) - do i=index+1,res_n-1 - if (res_tab(1,1,i)) then - done = res_tab(1,2,i) - else - res_tab(1,0,i) = done - res_tab(1,1,i) = done - res_tab(1,2,i) = done - endif - enddo - done = res_tab(1,0,index) - do i=index-1,0,-1 - if (res_tab(1,1,i)) then - done = res_tab(1,0,i) - else - res_tab(1,0,i) = done - res_tab(1,1,i) = done - res_tab(1,2,i) = done - endif - enddo - else - do i=0,res_n-1 - res_tab(1,0,i) = .true. - res_tab(1,1,i) = .true. - res_tab(1,2,i) = .true. - enddo - endif - -c Finally fill the last row with AND operation - do i=0,res_n-1 - do j=0,2 - res_tab(2,j,i) = (res_tab(0,j,i) .and. res_tab(1,j,i)) - enddo - enddo - - return - end - -c------------------------------------------------------------- - - subroutine construct_ranges(phi_n,phi_start,phi_end) -c Given the data in res_..., construct a table of -c min/max allowed angles - implicit none - -c Includes - include 'COMMON.GEO' - include 'COMMON.LOCMOVE' - -c Output arguments - integer phi_n - double precision phi_start(0:11),phi_end(0:11) - -c Local variables - logical done - integer index - - - if (res_n .eq. 0) then -c Any move is allowed - phi_n = 1 - phi_start(0) = -PI - phi_end(0) = PI - else - phi_n = 0 - index = 0 - done = .false. - do while (.not.done) -c Find start of range (01x) - done = .false. - do while (.not.done) - if (res_tab(2,0,index).or.(.not.res_tab(2,1,index))) then - index=index+1 - else - done = .true. - phi_start(phi_n) = res_ang(index) - endif - if (index .eq. res_n) done = .true. - enddo -c If a start was found (index < res_n), find the end of range (x10) -c It may not be found without wrapping around - if (index .lt. res_n) then - done = .false. - do while (.not.done) - if ((.not.res_tab(2,1,index)).or.res_tab(2,2,index)) then - index=index+1 - else - done = .true. - endif - if (index .eq. res_n) done = .true. - enddo - if (index .lt. res_n) then -c Found the end of the range - phi_end(phi_n) = res_ang(index) - phi_n=phi_n+1 - index=index+1 - if (index .eq. res_n) then - done = .true. - else - done = .false. - endif - else -c Need to wrap around - done = .true. - phi_end(phi_n) = flag - endif - endif - enddo -c Take care of the last one if need to wrap around - if (phi_end(phi_n) .eq. flag) then - index = 0 - do while ((.not.res_tab(2,1,index)).or.res_tab(2,2,index)) - index=index+1 - enddo - phi_end(phi_n) = res_ang(index) + 2.*PI - phi_n=phi_n+1 - endif - endif - - return - end - -c------------------------------------------------------------- - - subroutine fix_no_moves(phi) - implicit none - -c Includes - include 'COMMON.GEO' - include 'COMMON.LOCMOVE' - -c Output arguments - double precision phi - -c Local variables - integer index - double precision diff,temp - - -c Look for first 01x in gammas (there MUST be at least one) - diff = flag - index = 0 - do while (res_tab(1,0,index) .or. (.not.res_tab(1,1,index))) - index=index+1 - enddo - if (res_ang(index) .le. 0.D0) then ! Make sure it's from PHImax -c Try to increase PHImax - if (index .gt. 0) then - phi = res_ang(index-1) - diff = abs(res_ang(index) - res_ang(index-1)) - endif -c Look for last (corresponding) x10 - index = res_n - 1 - do while ((.not.res_tab(1,1,index)) .or. res_tab(1,2,index)) - index=index-1 - enddo - if (index .lt. res_n-1) then - temp = abs(res_ang(index) - res_ang(index+1)) - if (temp .lt. diff) then - phi = res_ang(index+1) - diff = temp - endif - endif - endif - -c If increasing PHImax didn't work, decreasing PHImin -c will (with one exception) -c Look for first x10 (there MUST be at least one) - index = 0 - do while ((.not.res_tab(1,1,index)) .or. res_tab(1,2,index)) - index=index+1 - enddo - if (res_ang(index) .lt. 0.D0) then ! Make sure it's from PHImin -c Try to decrease PHImin - if (index .lt. res_n-1) then - temp = abs(res_ang(index) - res_ang(index+1)) - if (res_ang(index+1) .le. 0.D0 .and. temp .lt. diff) then - phi = res_ang(index+1) - diff = temp - endif - endif -c Look for last (corresponding) 01x - index = res_n - 1 - do while (res_tab(1,0,index) .or. (.not.res_tab(1,1,index))) - index=index-1 - enddo - if (index .gt. 0) then - temp = abs(res_ang(index) - res_ang(index-1)) - if (res_ang(index-1) .ge. 0.D0 .and. temp .lt. diff) then - phi = res_ang(index-1) - diff = temp - endif - endif - endif - -c If it still didn't work, it must be PHImax == 0. or PHImin == PI - if (diff .eq. flag) then - index = 0 - if (res_tab(index,1,0) .or. (.not.res_tab(index,1,1)) .or. - + res_tab(index,1,2)) index = res_n - 1 -c This MUST work at this point - if (index .eq. 0) then - phi = res_ang(1) - else - phi = res_ang(index - 1) - endif - endif - - return - end - -c------------------------------------------------------------- - - integer function move_res(PHImin,PHImax,i_move) -c Moves residue i_move (in array c), leaving everything else fixed -c Starting geometry is not checked, it should be correct! -c R(,i_move) is the only residue that will move, but must have -c 1 < i_move < nres (i.e., cannot move ends) -c Whether any output is done is controlled by locmove_output -crc implicit none - -c Includes - implicit real*8 (a-h,o-z) - include 'DIMENSIONS' - include 'COMMON.CHAIN' - include 'COMMON.GEO' - include 'COMMON.LOCMOVE' - -c External functions - double precision ran_number - external ran_number - -c Input arguments - double precision PHImin,PHImax - integer i_move - -c RETURN VALUES: -c 0: move successfull -c 1: Dmin or Dmax had to be modified -c 2: move failed - check your input geometry - - -c Local variables - double precision X(0:2),Y(0:2),Z(0:2),Orig(0:2) - double precision P(0:2) - logical no_moves,done - integer index,i,j - double precision phi,temp,radius - double precision phi_start(0:11), phi_end(0:11) - integer phi_n - -c Set up the coordinate system - do i=0,2 - Orig(i)=0.5*(c(i+1,i_move-1)+c(i+1,i_move+1)) ! Position of origin - enddo - - do i=0,2 - Z(i)=c(i+1,i_move+1)-c(i+1,i_move-1) - enddo - temp=sqrt(Z(0)*Z(0)+Z(1)*Z(1)+Z(2)*Z(2)) - do i=0,2 - Z(i)=Z(i)/temp - enddo - - do i=0,2 - X(i)=c(i+1,i_move)-Orig(i) - enddo -c radius is the radius of the circle on which c(,i_move) can move - radius=sqrt(X(0)*X(0)+X(1)*X(1)+X(2)*X(2)) - do i=0,2 - X(i)=X(i)/radius - enddo - - Y(0)=Z(1)*X(2)-X(1)*Z(2) - Y(1)=X(0)*Z(2)-Z(0)*X(2) - Y(2)=Z(0)*X(1)-X(0)*Z(1) - -c Calculate min, max angles coming from dmin, dmax to c(,i_move-2) - if (i_move.gt.2) then - do i=0,2 - P(i)=c(i+1,i_move-2)-Orig(i) - enddo - call minmax_angles(P(0)*X(0)+P(1)*X(1)+P(2)*X(2), - + P(0)*Y(0)+P(1)*Y(1)+P(2)*Y(2), - + P(0)*Z(0)+P(1)*Z(1)+P(2)*Z(2), - + radius,a_n,a_ang,a_tab) - else - a_n=0 - endif - -c Calculate min, max angles coming from dmin, dmax to c(,i_move+2) - if (i_move.lt.nres-2) then - do i=0,2 - P(i)=c(i+1,i_move+2)-Orig(i) - enddo - call minmax_angles(P(0)*X(0)+P(1)*X(1)+P(2)*X(2), - + P(0)*Y(0)+P(1)*Y(1)+P(2)*Y(2), - + P(0)*Z(0)+P(1)*Z(1)+P(2)*Z(2), - + radius,b_n,b_ang,b_tab) - else - b_n=0 - endif - -c Construct the resulting table for alpha and beta - call construct_tab() - - if (locmove_output) then - print *,'ALPHAS & BETAS TABLE' - call output_tabs() - endif - -c Check that there is at least one possible move - no_moves = .true. - if (res_n .eq. 0) then - no_moves = .false. - else - index = 0 - do while ((index .lt. res_n) .and. no_moves) - if (res_tab(2,1,index)) no_moves = .false. - index=index+1 - enddo - endif - if (no_moves) then - if (locmove_output) print *,' *** Cannot move anywhere' - move_res=2 - return - endif - -c Transfer res_... into a_... - a_n = 0 - do i=0,res_n-1 - if ( (res_tab(2,0,i).neqv.res_tab(2,1,i)) .or. - + (res_tab(2,0,i).neqv.res_tab(2,2,i)) ) then - a_ang(a_n) = res_ang(i) - do j=0,2 - a_tab(j,a_n) = res_tab(2,j,i) - enddo - a_n=a_n+1 - endif - enddo - -c Check that the PHI's are within [0,PI] - if (PHImin .lt. 0. .or. abs(PHImin) .lt. small) PHImin = -flag - if (PHImin .gt. PI .or. abs(PHImin-PI) .lt. small) PHImin = PI - if (PHImax .gt. PI .or. abs(PHImax-PI) .lt. small) PHImax = flag - if (PHImax .lt. 0. .or. abs(PHImax) .lt. small) PHImax = 0. - if (PHImax .lt. PHImin) PHImax = PHImin -c Calculate min and max angles coming from PHImin and PHImax, -c and put them in b_... - call angles2tab(PHImin, PHImax, b_n, b_ang, b_tab) -c Construct the final table - call construct_tab() - - if (locmove_output) then - print *,'FINAL TABLE' - call output_tabs() - endif - -c Check that there is at least one possible move - no_moves = .true. - if (res_n .eq. 0) then - no_moves = .false. - else - index = 0 - do while ((index .lt. res_n) .and. no_moves) - if (res_tab(2,1,index)) no_moves = .false. - index=index+1 - enddo - endif - - if (no_moves) then -c Take care of the case where no solution exists... - call fix_no_moves(phi) - if (locmove_output) then - print *,' *** Had to modify PHImin or PHImax' - print *,'phi: ',phi*rad2deg - endif - move_res=1 - else -c ...or calculate the solution -c Construct phi_start/phi_end arrays - call construct_ranges(phi_n, phi_start, phi_end) -c Choose random angle phi in allowed range(s) - temp = 0. - do i=0,phi_n-1 - temp = temp + phi_end(i) - phi_start(i) - enddo - phi = ran_number(phi_start(0),phi_start(0)+temp) - index = 0 - done = .false. - do while (.not.done) - if (phi .lt. phi_end(index)) then - done = .true. - else - index=index+1 - endif - if (index .eq. phi_n) then - done = .true. - else if (.not.done) then - phi = phi + phi_start(index) - phi_end(index-1) - endif - enddo - if (index.eq.phi_n) phi=phi_end(phi_n-1) ! Fix numerical errors - if (phi .gt. PI) phi = phi-2.*PI - - if (locmove_output) then - print *,'ALLOWED RANGE(S)' - do i=0,phi_n-1 - print *,phi_start(i)*rad2deg,phi_end(i)*rad2deg - enddo - print *,'phi: ',phi*rad2deg - endif - move_res=0 - endif - -c Re-use radius as temp variable - temp=radius*cos(phi) - radius=radius*sin(phi) - do i=0,2 - c(i+1,i_move)=Orig(i)+temp*X(i)+radius*Y(i) - enddo - - return - end - -c------------------------------------------------------------- - - subroutine loc_test -crc implicit none - -c Includes - implicit real*8 (a-h,o-z) - include 'DIMENSIONS' - include 'COMMON.GEO' - include 'COMMON.LOCAL' - include 'COMMON.LOCMOVE' - -c External functions - integer move_res - external move_res - -c Local variables - integer i,j - integer phi_n - double precision phi_start(0:11),phi_end(0:11) - double precision phi - double precision R(0:2,0:5) - - locmove_output=.true. - -c call angles2tab(30.*deg2rad,70.*deg2rad,a_n,a_ang,a_tab) -c call angles2tab(80.*deg2rad,130.*deg2rad,b_n,b_ang,b_tab) -c call minmax_angles(0.D0,3.8D0,0.D0,3.8D0,b_n,b_ang,b_tab) -c call construct_tab -c call output_tabs - -c call construct_ranges(phi_n,phi_start,phi_end) -c do i=0,phi_n-1 -c print *,phi_start(i)*rad2deg,phi_end(i)*rad2deg -c enddo - -c call fix_no_moves(phi) -c print *,'NO MOVES FOUND, BEST PHI IS',phi*rad2deg - - R(0,0)=0.D0 - R(1,0)=0.D0 - R(2,0)=0.D0 - R(0,1)=0.D0 - R(1,1)=-cos(28.D0*deg2rad) - R(2,1)=-0.5D0-sin(28.D0*deg2rad) - R(0,2)=0.D0 - R(1,2)=0.D0 - R(2,2)=-0.5D0 - R(0,3)=cos(30.D0*deg2rad) - R(1,3)=0.D0 - R(2,3)=0.D0 - R(0,4)=0.D0 - R(1,4)=0.D0 - R(2,4)=0.5D0 - R(0,5)=0.D0 - R(1,5)=cos(26.D0*deg2rad) - R(2,5)=0.5D0+sin(26.D0*deg2rad) - do i=1,5 - do j=0,2 - R(j,i)=vbl*R(j,i) - enddo - enddo -c i=move_res(R(0,1),0.D0*deg2rad,180.D0*deg2rad) - imov=2 - i=move_res(0.D0*deg2rad,180.D0*deg2rad,imov) - print *,'RETURNED ',i - print *,(R(i,3)/vbl,i=0,2) - - return - end - -c-------------------------------------------------------------