{"id":312,"date":"2024-11-01T16:09:45","date_gmt":"2024-11-01T16:09:45","guid":{"rendered":"http:\/\/divergentfocus.co.uk\/?page_id=214"},"modified":"2024-11-01T16:09:45","modified_gmt":"2024-11-01T16:09:45","slug":"fortran-code","status":"publish","type":"page","link":"https:\/\/www.divergentfocus.co.uk\/?page_id=312","title":{"rendered":"Fortran Code"},"content":{"rendered":"\n

icorot.f<\/h2>\n\n\n\n


c icosahedron rotator program 0.1a
implicit real*8(a-h,o-z)
parameter(n=60)
double precision m,m2,xtemp,ytemp,ztemp,theta,a,b,p
double precision conv,switch,rad
dimension x(n),y(n),z(n)
open(12,file=’c60.xyz’)
open(13,file=’c60copy.xyz’)
open(14,file=’c60ico.xyz’)
switch=0.0d0
p=20.90515744515d0
10 format(1x,f24.20,1x,f24.20,1x,f24.20)
81 format(f24.20,1x,f24.20,1x,I10)
do i=1,n
read(12,10)x(i),y(i),z(i)
enddo
rad=3.1415926540d0\/180.0d0
do i=1,n
x(i)=x(i)*rad
y(i)=y(i)*rad
z(i)=z(i)*rad
enddo
do i=1,n
write(13,10)x(i),y(i),z(i)
enddo
do i=1,9999999
if (switch.eq.1) goto 91
p=p*rad+1.0E-14*rad
a=(y(1)\/rad)*sin(p)+(z(1)\/rad)*cos(p)
b=(y(2)\/rad)*sin(p)+(z(2)\/rad)*cos(p)
conv=abs(a-b)
p=p\/rad
c write(6,*)p,a,b,conv,i
write(6,81)p,conv,i
if (conv.lt.1.0E-15) switch=1
enddo
91 continue
p=p*rad
do i=1,n
ytemp=0.0d0
ztemp=0.0d0
ytemp=y(i)*cos(p)-z(i)*sin(p)
ztemp=y(i)*sin(p)+z(i)*cos(p)
y(i)=ytemp
z(i)=ztemp
enddo
write(13,*)’rotated’
do i=1,n
x(i)=x(i)\/rad
y(i)=y(i)\/rad
z(i)=z(i)\/rad
enddo
do i=1,n
write(13,10)x(i),y(i),z(i)
write(14,10)x(i),y(i),z(i)
enddo
stop
end<\/td><\/tr><\/tbody><\/table><\/figure>\n\n\n\n

pent.f<\/h2>\n\n\n\n
\n

 <\/p>\n\n \n
 <\/td>\n \n 
\nc      implicit real*8 (a-h,o-z)\n      character*4 itex,itexh,itexb,ibas,ibash,istg\n      integer p, s\n      double precision aconst,r,tconst,zdf,test,pol,rr,dota,dotb\n      double precision ratph, rstart, rstep, rpent, rhex,switch\n      double precision conv\n      parameter (n=5)\n      parameter (s=200)\n      dimension x(12),y(12),z(12),xc(20),yc(20),zc(20)\n      dimension xic(30),yic(30),zic(30),xp(12),yp(12),op(12)\n      dimension xb(12),yb(12),zb(12),xbh(12),ybh(12),zbh(12)\n      dimension xh(20),yh(20),zh(20),xi(60),yi(60),zi(60)\n      Dimension F(3), alp(3,3), xyz(n,3), b(3), d(3), t(3,3)\n      dimension af(3), at(3,3), stmew(n,3), cor(n,n,3), old(3)\n      dimension c(3), e(3), g(3), stold(n,3), pent(5)\n      dimension vx(n), vy(n), vz(n),stdif(n,3)\n      data xp\/0.5d0,-0.5d0,0.0d0,5.0d0,5.0d0,0.0d0,-5.0d0,0.0d0,\n     10.5d0,0.0d0,-5.0d0,-0.5d0\/\n      data yp\/0.0d0,0.0d0,5.0d0,0.5d0,-0.5d0,-5.0d0,0.5d0,5.0d0,\n     10.0d0,-5.0d0,-0.5d0,0.0d0\/\n      data op\/5.0d0,5.0d0,0.5d0,0.0d0,0.0d0,0.5d0,0.0d0,-0.5d0,\n     1-5.0d0,-0.5d0,0.0d0,-5.0d0\/\n      data itex\/'C   '\/\n      data itexh\/'H   '\/\n      data itexb\/'B   '\/\n      data ibas\/'P   '\/\n      data ibash\/'S   '\/\n      data istg\/' 12C'\/\n      data conv\/0.5291771d0\/\n      zdf=0.1d0\n      aconst=0.25d0\n      m=0.0d0\n      switch=0.0d0\n      conv=0.0d0\n      open(14,file='cplay.xyz')\n      open(15,file='q.in')\n      open(13,file='c60.xyz')\n      open(7,file='f.dat')\n      open(8,file='data.out')\n      open(10,file='qrun.out')\n      open(12,file='debug.dat')\n      open(16,file='c60.dat')\n      open(17,file='spread.out')\n\n      read(15,*)ratph,rstart,rstep\n      write(6,*)'enter pent numbers'\n      read(5,*)pent\n\nc     The q-loop!\n      do q=1,1\n      do i=1,12\n         x(i)=xp(i)\n         y(i)=yp(i)\n         z(i)=op(i)\n      enddo\n\nc     chris here\n      rewind 16\n      rewind 7\n\n      if (q.ne.1) rstart=rstart+rstep\n      rpent=rstart\n      rhex=ratph*rpent\n      write(10,*)'ratph',ratph,'rstart',rstart\n      write(10,*)'rpent',rpent,'rhex',rhex,'rstep',rstep\n      read(16,*)rbb,rbh\n      read(16,*)rcc,rch\n      write(10,*)'rbb',rbb,'rbh',rbh\n      write(10,*)'rcc',rcc,'rch',rch\nc     chris ends here\n\n      rbb=rbb\/conv\n      rbh=rbh\/conv\n      rcc=rcc\/conv\n      rch=rch\/conv\n      rpent=rpent\/conv\n      rhex=rhex\/conv\n      char=6.0d0\n      charh=1.0d0\n      charb=5.0d0\n      phi=0.5d0*(1.0d0+dsqrt(5.0d0))\n      i=0\n      j=0\n      k=0\n      do 1 i=1,12\n      if(x(i).gt.2.0d0)x(i)=phi\/2.0d0\n      if(x(i).lt.-2.0d0)x(i)=-phi\/2.0d0\n      if(y(i).gt.2.0d0)y(i)=phi\/2.0d0\n      if(y(i).lt.-2.0d0)y(i)=-phi\/2.0d0\n      if(z(i).gt.2.0d0)z(i)=phi\/2.0d0\n      if(z(i).lt.-2.0d0)z(i)=-phi\/2.0d0\n1     continue\nc     read the faces of the icosahedron\n      do 2 ic=1,20\n      read(16,*)i, j, k\n      xc(ic)=(x(i)+x(j)+x(k))\/phi\n      yc(ic)=(y(i)+y(j)+y(k))\/phi\n      zc(ic)=(z(i)+z(j)+z(k))\/phi\n2     continue\n\nc      write(10,*)'x(1) at 2',x(1)\n\n      write(6,89)\n89    format(1X,' WARNING!! INPUT ANGSTROM OUTPUT BOHR!!!!! ')\n      write(6,90)\n90    format(1X,' ICOSAHEDRON OF UNIT SIDE ')\n      do 3 i=1,9\n3     write(6,100)itex,i,ibas,x(i),y(i),z(i),char\n      do 33 i=10,12\n33    write(6,101)itex,i,ibas,x(i),y(i),z(i),char\n100   format(A1,I1,2X,2X,A2,3X,4F12.9)\n101   format(A1,I2,1X,2X,A2,3X,4F12.9)\n      write(6,80)rbb\n80    format(1X,' B12 POSITIONS FOR RBB= ',1X,F10.6)\nc     scale for boron cage of B12H12\n\nc     ze problem is ze ere\n\n      do 81 i=1,12\n      xb(i)=x(i)*rbb\n      yb(i)=y(i)*rbb\n      zb(i)=z(i)*rbb\n      if(i.le.9)then\n      write(6,100)itexb,i,ibas,xb(i),yb(i),zb(i),charb\n      else\n      write(6,101)itexb,i,ibas,xb(i),yb(i),zb(i),charb\n      endif\n81    continue\nc     scale for H cage\n      rad=dsqrt(xb(1)**2.0d0+yb(1)**2.0d0+zb(1)**2.0d0)\n\nc      JIZZ!! taking off 'd' what will it do??\n\n      write(6,82)rbb,rbh\n82    format(1X,' H12 POSITIONS FOR RBB,RBH= ',1X,F10.6,1X,F10.6)\n      do 83 i=1,12\n      xbh(i)=x(i)*rbb*(rbh+rad)\/rad\n      ybh(i)=y(i)*rbb*(rbh+rad)\/rad\n      zbh(i)=z(i)*rbb*(rbh+rad)\/rad\n      if(i.le.9)then\n      write(6,100)itexh,i,ibash,xbh(i),ybh(i),zbh(i),charh\n      else\n      write(6,101)itexh,i,ibash,xbh(i),ybh(i),zbh(i),charh\n      endif\n83    continue\n      write(6,91)\n91    format(1X,' DODECAHEDRON OF UNIT SIDE ')\n      do 4 i=1,9\n4     write(6,100)itex,i,ibas,xc(i),yc(i),zc(i),char\n      do 34 i=10,20\n34    write(6,101)itex,i,ibas,xc(i),yc(i),zc(i),char\n      write(6,92)rcc\n92    format(1X,' C20 WITH SIDE RCC=',1X,F10.6)\nc     scale dodecahedron\n      do 5 i=1,20\n      xc(i)=rcc*xc(i)\n      yc(i)=rcc*yc(i)\n      zc(i)=rcc*zc(i)\n      if(i.le.9)then\n      write(6,100)itex,i,ibas,xc(i),yc(i),zc(i),char\n      else\n      write(6,101)itex,i,ibas,xc(i),yc(i),zc(i),char\n      endif\n5     continue\n\n      rad=dsqrt(xc(1)**2.0d0+yc(1)**2.0d0+zc(1)**2.0d0)\n      write(6,93)rcc,rch\n93    format(1X,' H20 POSITIONS FOR RCC,RCH =',1X,F10.6,1X,F10.6)\n      do 6 i=1,20\n      xh(i)=xc(i)*(rch+rad)\/rad\n      yh(i)=yc(i)*(rch+rad)\/rad\n      zh(i)=zc(i)*(rch+rad)\/rad\n      if(i.le.9)then\n      write(6,100)itexh,i,ibash,xh(i),yh(i),zh(i),charh\n      else\n      write(6,101)itexh,i,ibash,xh(i),yh(i),zh(i),charh\n      endif\n6     continue\n      write(12,94)rpent,rhex\n94    format(1X,' C60 POSITIONS FOR RPENT,RHEX = ',2(1X,F10.6))\n      alpha=rpent\/(2.0D0*rpent+rhex)\n      rbig=rhex+2.0d0*rpent\nc     scale master icosahedron\n      do 7 i=1,12\n      x(i)=x(i)*rbig\n      y(i)=y(i)*rbig\n      z(i)=z(i)*rbig\n7     continue\nc     read edges of icosahedron\n      do 8 i=1,30\n      read(16,*)j,k\n      xi(i)=x(j)+alpha*(x(k)-x(j))\n      yi(i)=y(j)+alpha*(y(k)-y(j))\n      zi(i)=z(j)+alpha*(z(k)-z(j))\n      xi(i+30)=x(k)+alpha*(x(j)-x(k))\n      yi(i+30)=y(k)+alpha*(y(j)-y(k))\n      zi(i+30)=z(k)+alpha*(z(j)-z(k))\n      xic(i)=x(j)+x(k)\n      yic(i)=y(j)+y(k)\n      zic(i)=z(j)+z(k)\n      ric=xic(i)*xic(i)+yic(i)*yic(i)+zic(i)*zic(i)\n      xic(i)=xic(i)\/ric\n      yic(i)=yic(i)\/ric\n      zic(i)=zic(i)\/ric\n8     continue\n\nc     altered...\n123   format(4x,d20.16,1x,d20.16,1x,d20.16)\n      do 888 i=1,9\nc 888   write(12,100)itex,i,ibas,xi(i),yi(i),zi(i),char\n888   write(14,123)xi(i),yi(i),zi(i)\n      do 889 i=10,60\nc 889   write(12,101)itex,i,ibas,xi(i),yi(i),zi(i),char\n889   write(14,123)xi(i),yi(i),zi(i)\n\nc     altered... \n      write(10,*)pent\n      vx(1)=xi(pent(1))\n      vy(1)=yi(pent(1))\n      vz(1)=zi(pent(1))\n      vx(2)=xi(pent(2))\n      vy(2)=yi(pent(2))\n      vz(2)=zi(pent(2))\n      vx(3)=xi(pent(3))\n      vy(3)=yi(pent(3))\n      vz(3)=zi(pent(3))\n      vx(4)=xi(pent(4))\n      vy(4)=yi(pent(4))\n      vz(4)=zi(pent(4))\n      vx(5)=xi(pent(5))\n      vy(5)=yi(pent(5)) \n      vz(5)=zi(pent(5))\n\n      zero=0.0D0\n      izero=0\n      do 977 i=1,60\n      xi(i)=xi(i)*conv\n      yi(i)=yi(i)*conv\n977   zi(i)=zi(i)*conv\n      do 988 i=1,60\n988   write(6,111)istg,izero,izero,xi(i),yi(i),zi(i),izero,izero,zero\n111   format(1X,A4,2I5,3F10.7,2I5,F10.4)\n      write(6,95)\n95    format(1X,' ICOSIDODECAHEDRON OF UNIT RADIUS ')\n      do 444 i=1,30\n444   write(6,445)i,xic(i),yic(i),zic(i)\n445   format(1x,i3,3f16.10)\nc     polar coordinates of c60 in this orientation\n      write(6,2222)\n2222  format(1X,' POLAR COORDINATES OF C60 ')\n      pi=3.1415926536d0\n      do 989 i=1,60\n      r=sqrt(xi(i)*xi(i)+yi(i)*yi(i)+zi(i)*zi(i))\n      xi(i)=xi(i)\/r\n      yi(i)=yi(i)\/r\n      zi(i)=zi(i)\/r\n      temp=zi(i)\n      if(temp.gt.1.0d0)temp=1.0d0\n      if(temp.lt.-1.0d0)temp=-1.0d0\n      tang=acos(temp)\n      s1=sin(tang)\n      test=s1*s1\n      if(test.lt.1.0d-6)goto 2727\n      temp=xi(i)\/s1\n      if(temp.gt.1.0d0)temp=1.0d0\n      if(temp.lt.-1.0d0)temp=-1.0d0\n      ping=acos(temp)\n      goto 2828\n2727  ping=0.0d0\n2828  continue\n      if(yi(i).lt.0.0d0)ping=2.0d0*pi-ping\n      tang=tang*180.0d0\/pi\n      ping=ping*180.0d0\/pi\n      write(6,1111)i,xi(i),yi(i),zi(i),tang,ping\n1111  format(i3,5f14.4)\n989   continue\n\nc     chris's program starts here\nc     *****************************************************\n      do i=1,3\n        read(7,*)f(i)\n        write(6,*)'reading f from f.dat',f(i)\n      enddo\n 10   format(4x,d20.16,1x,d20.16,1x,d20.16)\n      rr=SQRT((vx(1)**2)+(vy(1)**2)+(vz(1)**2))\n      do i=1,n\n      xyz(i,1)=vx(i)\n      xyz(i,2)=vy(i)\n      xyz(i,3)=vz(i)\n      enddo\nc     *****************************************************\nc     set up alpha\n      do i=1,3\n        do j=1,3\n          if (i.eq.j) then\n            alp(i,j)=aconst\n          else\n            alp(i,j)=0.0d0\n          endif\n        enddo\n      enddo\nc     ******************************************************\n      zdff=1.0d0\nc     h loop numbet of total cycles set by 's'\n      do h=1,10\n      if (h.ne.1) zdff=zdf\n      if(switch.eq.1) goto 67\n      do l=1,n\n      do m=1,n\n         if (m.eq.l) then\n           goto 98\n         else\n         endif\nc     ******************************************************\nc     d(3) displacement vector between l and m\n      do i=1,3\n        d(i) = 0.0d0\n        d(i) = xyz(l,i)-xyz(m,i)\n      enddo\nc     r distance between i and j\n      r=SQRT ((d(1)**2)+(d(2)**2)+(d(3)**2))\nc     af multiplying alp(3,3) by F(3) to get af(3)\n      do i=1,3\n        af(i)=0.0d0\n        do j=1,3\n          af(i)=af(i)+(alp(i,j)*(f(j)))\n        enddo\n      enddo\nc     t construction\n      tconst=3.0d0\/(r**5)\n      do i=1,3\n        do j=1,3\n          if (i.eq.j) then\n            t(i,j)=tconst*((D(i)**2)-((1.0d0\/3.0d0)*(r**2)))\n          else\n            t(i,j)=tconst*(d(i)*(d(j)))\n          end if\n        enddo\n      enddo\nc     at construction\n      do i=1,3\n        do j=1,3\n          at(i,j)=0.0d0\n          do k=1,3\n            at(i,j)=at(i,j)+alp(i,k)*t(k,j)\n          enddo\n        enddo\n      enddo\nc     *****************************************************\n      do i=1,3\n        if (h.eq.1) then\n          stmew(l,i)=af(i)\n          old(i)=stmew(l,i)\n        else\n         old(i)=stmew(l,i)\n        endif\n      enddo\n      do i=1,3\n          cor(l,m,i)=0.0d0\n        do j=1,3\n          cor(l,m,i)=cor(l,m,i)+(at(i,j)*(old(j)))\n        enddo\n      enddo\n  98  continue\n      enddo\nc      write(10,*)'r',r\n      enddo\nc     ************* end l and m loops *********************\n      do l=1,n\n        do i=1,3\n          stmew(l,i)=af(i)\n          if (h.eq.1) stold(l,i)=af(i)\n          do j=1,n\n            if (j.eq.l) then\n              write(10,*)h,conv\n            else\n              stmew(l,i)=stmew(l,i)+cor(j,l,i)\n              write(6,*)stmew(l,i)\n            endif\n           enddo\n        enddo\n      enddo\n      do l=1,n\n        do i=1,3\n          stmew(l,i)=zdff*stmew(l,i)+(1.0d0-zdff)*stold(l,i)\n          stdif(l,i)=abs(stmew(l,i)-stold(l,i))\n          dota=dota+stdif(l,i)*stdif(l,i)\n          dotb=dotb+stmew(l,i)*stmew(l,i)\n          conv=dota\/dotb\n          if (conv.lt.0.001d0) switch=1\n          stold(l,i)=stmew(l,i)\n        enddo\n      enddo\n      enddo\n   67 continue\nc     ************ end of h loop **************************\n  110 do l=1,n\nc        write(12,*)'particle',l\nc        do i=1,3\nc           write(12,*)stmew(l,i)\nc        enddo\n      enddo\nc     finding polarisability\n      pol=0.0d0\n      do l=1,n\n      pol=pol+stmew(l,3)\n      enddo\n      write(12,*)'r-displacement from origin',rr\n      write(12,*)'sum dipoles ',pol\n      write(12,*)'field ', f(3)\n      pol=pol\/f(3)\n      write(12,*)'polarisability',pol\n      write(10,*)'polarisability',pol\n      write(10,*)'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'\n      write(6,*)'polarisability',pol\n      rpent=rpent*conv\n      rhex=rhex*conv\n      write(17,45)rpent,rhex,rr,pol\n  45  format(f8.4,3x,f14.8,3x,d24.8,2x,d23.14)\n  20  format(2x,f23.20,1x,f23.20,1x,f23.20)\n      enddo\n  99  stop\n      end\n \n\n  <\/pre>\n      <\/td>\n  <\/tr>\n<\/table>\n\n\n\n
vsh.f<\/h2>\n\n\n\n

\nc     VSH 0.9d\n      implicit real*8 (a-h,o-z)\n      parameter (n=60)\n      double precision soconst, v4k, v5k, r, kk, dotp,v6k,v7k,v8k\n      double precision v9k,v100k,v110k,k0,k1,k5,k7,k9,k11,k13\n      double precision v130k,v3k,k3,v150k,k15,k17\n      double precision a1,a2,a3,a4,a5,a6,a7,a8,a9\n      dimension xyz(n,3),x(n),y(n),z(n),so(n,3),v10(n,3),v20(n,3)\n      dimension v30(n,3), v40(n,3), v50(n,3), vnorm(n,3), c60(n,3)\n      dimension vorth(n,3),xi(n),yi(n),zi(n),v60(n,3),v70(n,3),v80(n,3)\n      dimension v90(n,3),v100(n,3),v110(n,3),cage(n,3),v130(n,3)\n      dimension sum(3),v150(n,3),v170(n,3)\n      open(14,file='test.out')\nc      open(13,file='c60ico.xyz')\n      open(13,file='c60.xyz')\n      open(12,file='vsh.out')\n      open(11,file='vsh.in')\n      open(10,file='check.out')\n      open(9,file='vshc60.out')\nc     **********************************************************\nc     reading in co-ordinates of c60 from c60.xyz\n      do i=1,n\n      read(13,10)x(i), y(i), z(i)\n      enddo\nc     reading in dipoles from vsh.in\n      a1=1.0d0\n      write(14,28)a1,a1,a1\n      do i=1,n\n      read(11,28)xi(i),yi(i),zi(i)\n      c60(i,1)=xi(i)\n      c60(i,2)=yi(i)\n      c60(i,3)=zi(i)\n      enddo\n   10 format(1x,d24.20,1x,d24.20,1x,d24.20)\n   40 format(2x,f23.20,1x,f23.20,1x,f23.20)\n   50 format(2x,f23.20,2x,f23.20,2x,f23.20)\n   25 format(4x,d20.12,1x,d20.12,1x,d20.12)\n   28 format(1x,d24.16,1x,d24.16,1x,d24.16)\n      r=SQRT ((x(1)**2)+(y(1)**2)+(z(1)**2))\n      write(6,*)'r',r\n      r=SQRT ((x(2)**2)+(y(2)**2)+(z(2)**2))\n      write(6,*)'r',r\n      r=SQRT ((x(3)**3)+(y(3)**2)+(z(3)**2))\n      write(6,*)'r',r\n      do i=1,n\n         x(i)=x(i)\/r\n         y(i)=y(i)\/r\n         z(i)=z(i)\/r\n      enddo\n      do i=1,10\n      write(6,*)'z',z(i)\n      enddo\nc     **********************************************************\n      write(12,*)'c60 dipole field'\n      call output(c60,n)\n\nc     So (0,0,1)\n      write(12,*)'So'\nc      soconst=1.0d0\/(60.0d0**0.5d0)\n      do i=1,n\n         so(i,1)=0.0d0\n         so(i,2)=0.0d0\n         so(i,3)=1.0d0\n      enddo\n      call normal(vnorm,n,so)\n      call vector(n,so,vnorm)\n      call output(so,n)\n\nc      orthog c60 to so here!\n      call orthog(vnorm,n,c60,so)\n      call vector(n,c60,vnorm)\n\n      call checkn(so,n)\n  \n      call dotprod(k0,n,so,c60)\n      write(9,*)'dotp so c60',k0\n      call cv(so,n)\n      write(12,*)'c60 dipole field after so'\n      call output(c60,n)\n      call cv(c60,n)\n      call sumit(n,c60,sumo)\nc     **********************************************************\nc     V10 (-xz, -yz, 1-zz)\n      write(12,*)'V10'\n      do i=1,n\n        v10(i,1)=(-1.0d0*(x(i)*z(i)))\n        v10(i,2)=(-1.0d0*(y(i)*z(i)))\n        v10(i,3)=1.0d0-((z(i)*z(i)))\n      enddo\n   30 format(f23.20)\n      write(10,*)'V10 raw'\n      call checko(n,v10,so)\n      write(10,*)'v10 end raw'\n      call output(v10,n)\n  \n      call normal(vnorm,n,v10)\n      call vector(n,v10,vnorm)\n      call orthog(vnorm,n,v10,so)\n      call vector(n,v10,vnorm)\n      call normal(vnorm,n,v10)\n      call vector(n,v10,vnorm)\n\nc     v10 orth to c60 here!\n      call orthog(vnorm,n,c60,v10)\n      call vector(n,c60,vnorm)\n\n      call checkn(v10,n)\n      call checko(n,v10,so)\n      call dotprod(k1,n,v10,c60)\n      write(9,*)'dotp v10 c60',k1\nc      call cv(v10,n)\n      write(12,*)'c60 dipole field after v10'\n      call output(c60,n)\n      call cv(c60,n)\n      call sumit(n,c60,sumo)\nc     ***********************************************************\nc      V20 3(-xz2, -yz2, z(1-z2)\n      do i=1,n\n         v20(i,1)=3.0d0*(-x(i)*(z(i)**2))\n         v20(i,2)=3.0d0*(-y(i)*(z(i)**2))\n         v20(i,3)=3.0d0*(z(i)*(1.0d0-(z(i)**2)))\n      enddo\n      write(12,*)'V20'\nc      call output(v20,n)\n      write(10,*)'V20 raw'\n      call checko(n,v20,so)\n      call checko(n,v20,v10)\n      write(10,*)'end raw'\n    \n      call normal(vnorm,n,v20) \n      call vector(n,v20,vnorm)\n\n      call checkn(v20,n)\n      call checko(n,v20,so)\n      call checko(n,v20,v10)\n\n      call dotprod(dotp,n,v20,c60)\n      write(9,*)'dotp v20 c60',dotp\nc     ***********************************************************\nc     V30\n      do i=1,n\n         v3k=(15.0d0*(z(i)**2.0d0))-3.0d0\n         v30(i,1)=-0.5d0*x(i)*z(i)*v3k\n         v30(i,2)=-0.5d0*y(i)*z(i)*v3k\n         v30(i,3)=0.5d0*(1.0d0-(z(i)**2.0d0))*v3k\n      enddo\n      write(12,*)'V30'\n      call output(v30,n)\n      write(10,*)'V30 raw'\n      call dotprod(dotp,n,v30,c60)\n      write(6,*)dotp\n      call checko(n,v30,so)\n      call checko(n,v30,v10)\n      call checko(n,v30,v20)\n      write(10,*)'v30 end raw'\n\n      call orthog(vnorm,n,v30,so)\n      call vector(n,v30,vnorm)\n      call normal(vnorm,n,v30) \n      call vector(n,v30,vnorm)\n      call orthog(vnorm,n,v30,v10)\n      call vector(n,v30,vnorm)\n      call normal(vnorm,n,v30) \n      call vector(n,v30,vnorm)\n\nc     make orthog v30 to c60 here! \n      call orthog(vorth,n,c60,v30)\n      call vector(n,c60,vorth)\n\n      call checko(n,v30,so)\n      call checko(n,v30,v10)\n      call checko(n,v30,v20)\n      call dotprod(k3,n,v30,c60)\n      write(9,*)'dotp v30 c60',k3\nc      call cv(v30,n)\n      write(12,*)'c60 dipole field after v30'\n      call output(c60,n)\n      call cv(c60,n)\n      call sumit(n,c60,sumo)\nc     ***********************************************************\nc     V40\n      v4k=0.0d0\n      do i=1,n\n         v4k=((35.0d0*(z(i)**3.0d0))-(15.0d0*z(i)))\n         v40(i,1)=-0.5d0*x(i)*z(i)*v4k\n         v40(i,2)=-0.5d0*y(i)*z(i)*v4k\n         v40(i,3)=0.5d0*(1.0d0-(z(i)*z(i)))*v4k\n      enddo\nc      write(12,*)'V40'\nc      call output(v40,n)\n      write(10,*)'v40 raw'\n      call checko(n,v40,so)\n      call checko(n,v40,v10)\n      call checko(n,v40,v20)\n      call checko(n,v40,v30)\n      write(10,*)'v40 raw end'\n  \n      call normal(vnorm,n,v40) \n      call vector(n,v40,vnorm)\nc      call orthog(vnorm,n,c60,v40)\nc      call vector(n,c60,vnorm)\nc      call checkn(c60,n)\nc      call checko(n,c60,v40)\n\nc      call orthog(vnorm,n,v40,v20)\nc      call vector(n,v40,vnorm)\nc      call checkn(v40,n)\n      call checko(n,v40,so)\n      call checko(n,v40,v10)\n      call checko(n,v40,v20)\n      call checko(n,v40,v30)\n\n      call dotprod(dotp,n,v40,c60)\n      write(9,*)'dotp v40 c60',dotp\nc     ***********************************************************\nc     V50\n      do i=1,n\n         v5k=(315.0d0*(z(i)**4.0d0))-(210.0d0*(z(i)**2.0d0))+15.0d0\n         v50(i,1)=-0.125d0*x(i)*z(i)*v5k\n         v50(i,2)=-0.125d0*y(i)*z(i)*v5k\n         v50(i,3)=0.125d0*(1.0d0-(z(i)*z(i)))*v5k\n      enddo\n      write(10,*)'V50raw'\n      write(12,*)'V50'\n      call output(v50,n)\n      call checko(n,v50,so)\n      call checko(n,v50,v10)\n      call checko(n,v50,v20)\n      call checko(n,v50,v30)\n      call checko(n,v50,v40)\n      write(10,*)'v50 end raw'\n\n      call orthog(vnorm,n,v50,so)\n      call vector(n,v50,vnorm)\n      call normal(vnorm,n,v50) \n      call vector(n,v50,vnorm)\n      call orthog(vnorm,n,v50,v10)\n      call vector(n,v50,vnorm)\n      call normal(vnorm,n,v50) \n      call vector(n,v50,vnorm)\n      call orthog(vnorm,n,v50,v30)\n      call vector(n,v50,vnorm)\n      call normal(vnorm,n,v50) \n      call vector(n,v50,vnorm)\n\nc     make orthog v50 to c60 here! \n      call orthog(vorth,n,c60,v50)\n      call vector(n,c60,vorth)\n\n      call checko(n,v50,so)\n      call checko(n,v50,v10)\n      call checko(n,v50,v20)\n      call checko(n,v50,v30)\n      call checko(n,v50,v40)\n      call dotprod(k5,n,v50,c60)\n      write(9,*)'dotp v50 c60',k5\nc      call cv(v50,n)\n      write(12,*)'c60 dipole field after v50'\n      call output(c60,n)\nc      c60(24,3)=0.0d0\nc      c60(25,3)=0.0d0\nc      c60(29,3)=0.0d0\nc      c60(32,3)=0.0d0\nc      c60(35,3)=0.0d0\nc      c60(36,3)=0.0d0\nc      c60(37,3)=0.0d0\nc      c60(56,3)=0.0d0\nc      write(12,*)'c60 dipole field after v50'\nc      call output(c60,n)\n      call cv(c60,n)\n      call sumit(n,c60,sumo)\nc     *******************************************************\nc     v60\n      write(12,*)'V60'\n      do i=1,n\n         v6k=((693.0d0*(z(i)**5.0d0))-(630.0d0*(z(i)**3.0d0)))\n         v6k=v6k+(105.0d0*z(i))\n         v60(i,1)=(-0.125d0*x(i)*z(i))*v6k\n         v60(i,2)=(-0.125d0*y(i)*z(i))*v6k\n         v60(i,3)=(0.125d0*(1.0d0-(z(i)**2)))*v6k\n      enddo\nc      call output(v60,n)\n      write(10,*)'V60 raw'\n      call checko(n,v60,so)\n      call checko(n,v60,v10)\n      call checko(n,v60,v20)\n      call checko(n,v60,v30)\n      call checko(n,v60,v40)\n      call checko(n,v60,v50)\n      write(10,*)'v60 end raw'\n\n      call orthog(vnorm,n,v60,v20)\n      call vector(n,v60,vnorm)\n      call checkn(v60,n)\n      call orthog(vnorm,n,v60,v40)\n      call vector(n,v60,vnorm)\n      call checkn(v60,n)\n\n\n      call checko(n,v60,so)\n      call checko(n,v60,v10)\n      call checko(n,v60,v20)\n      call checko(n,v60,v30)\n      call checko(n,v60,v40)\n      call checko(n,v60,v50)\n\n      call dotprod(dotp,n,v60,c60)\n      write(9,*)'dotp v60 c60',dotp\nc     ***********************************************************\nc     v70\n      write(12,*)'V70'\n      do i=1,n\n         v7k=(3003.0d0*(z(i)**6.0d0))-(3465.0d0*(z(i)**4.0d0))\n         v7k=v7k+(945.0d0*(z(i)**2.0d0))-35.0d0\n         v70(i,1)=-0.0625d0*x(i)*z(i)*v7k\n         v70(i,2)=-0.0625d0*y(i)*z(i)*v7k\n         v70(i,3)=(0.0625d0*(1.0d0-z(i)**2))*v7k\n      enddo\n      call output(v70,n)\n      write(10,*)'V70 raw'\n      call checko(n,v70,so)\n      call checko(n,v70,v10)\n      call checko(n,v70,v20)\n      call checko(n,v70,v30)\n      call checko(n,v70,v40)\n      call checko(n,v70,v50)\n      call checko(n,v70,v60)\n      write(10,*)'v70 end raw'\n\n      call orthog(vnorm,n,v70,so)\n      call vector(n,v70,vnorm)\n      call normal(vnorm,n,v70) \n      call vector(n,v70,vnorm)\n      call orthog(vnorm,n,v70,v10)\n      call vector(n,v70,vnorm)\n      call normal(vnorm,n,v70) \n      call vector(n,v70,vnorm)\n      call orthog(vnorm,n,v70,v30)\n      call vector(n,v70,vnorm)\n      call normal(vnorm,n,v70) \n      call vector(n,v70,vnorm)\n      call orthog(vnorm,n,v70,v50)\n      call vector(n,v70,vnorm)\n      call normal(vnorm,n,v70) \n      call vector(n,v70,vnorm)\n\nc     set to orthog c60 to v70\n      call orthog(vnorm,n,c60,v70)\n      call vector(n,c60,vnorm)\n \n      call checko(n,v70,so)\n      call checko(n,v70,v10)\n      call checko(n,v70,v20)\n      call checko(n,v70,v30)\n      call checko(n,v70,v40)\n      call checko(n,v70,v50)\n      call checko(n,v70,v60)\n      call dotprod(k7,n,v70,c60)\n      write(9,*)'dotp v70 c60',k7\nc     call cv(v70,n)\n      write(12,*)'c60 dipole field after v70'\n      call output(c60,n)\n      call cv(c60,n)\n      call sumit(n,c60,sumo)\nc     ***********************************************************\nc     v80\n      write(12,*)'V80'\n      do i=1,n\n         v8k=(6435*(z(i)**7))-(9009*(z(i)**5))+(3465*(z(i)**3))\n         v8k=v8k-(315*z(i))\n         v80(i,1)=-0.0625d0*x(i)*z(i)*v8k\n         v80(i,2)=-0.0625d0*y(i)*z(i)*v8k\n         v80(i,3)=(0.0625d0*(1.0d0-z(i)**2))*v8k\n      enddo\n      call output(v80,n)\n      write(10,*)'V80 raw'\n      call checko(n,v80,so)\n      call checko(n,v80,v10)\n      call checko(n,v80,v20)\n      call checko(n,v80,v30)\n      call checko(n,v80,v40)\n      call checko(n,v80,v50)\n      call checko(n,v80,v60)\n      call checko(n,v80,v70)\n      write(10,*)'v80 end raw'\n\n      call orthog(vnorm,n,v80,v20)\n      call vector(n,v80,vnorm)\n      call checkn(v80,n)\n      call orthog(vnorm,n,v80,v40)\n      call vector(n,v80,vnorm)\n      call checkn(v80,n)\n      call orthog(vnorm,n,v80,v60)\n      call vector(n,v80,vnorm)\n      call checkn(v80,n)\n\n      call checko(n,v80,so)\n      call checko(n,v80,v10)\n      call checko(n,v80,v20)\n      call checko(n,v80,v30)\n      call checko(n,v80,v40)\n      call checko(n,v80,v50)\n      call checko(n,v80,v60)\n      call checko(n,v80,v70)\n\n      call dotprod(dotp,n,v80,c60)\n      write(9,*)'dotp v80 c60',dotp\nc     ************************************************************\nc     v90\n      write(12,*)'v90'\nc      x(1)=4.0d0\nc      y(1)=3.0d0\nc      z(1)=2.0d0\n      do i=1,n\n         v9k=(109395.0d0*(z(i)**8.0d0))-(180180.0d0*(z(i)**6.0d0))\n         v9k=v9k+(90090.0d0*(z(i)**4.0d0))\n         v9k=v9k-(13860.0d0*(z(i)**2.0d0))+315.0d0\n         v90(i,1)=-(1.0d0\/128.0d0)*x(i)*z(i)*v9k\n         v90(i,2)=-(1.0d0\/128.0d0)*y(i)*z(i)*v9k\n         v90(i,3)=(1.0d0\/128.0d0)*(1.0d0-z(i)**2)*v9k\n      enddo\nc      write(6,*)v90(1,1),v90(1,2),v90(1,3)\n      call output(v90,n)\n      write(10,*)'v90 raw'\n      call checko(n,v90,so)\n      call checko(n,v90,v10)\n      call checko(n,v90,v20)\n      call checko(n,v90,v30)\n      call checko(n,v90,v40)\n      call checko(n,v90,v50)\n      call checko(n,v90,v60)\n      call checko(n,v90,v70)\n      call checko(n,v90,v80)\n      write(10,*)'v90 end raw'\n\n      call orthog(vnorm,n,v90,so)\n      call vector(n,v90,vnorm)\n      call normal(vnorm,n,v90) \n      call vector(n,v90,vnorm)\n      call orthog(vnorm,n,v90,v10)\n      call vector(n,v90,vnorm)\n      call normal(vnorm,n,v90) \n      call vector(n,v90,vnorm)\n      call orthog(vnorm,n,v90,v30)\n      call vector(n,v90,vnorm)\n      call normal(vnorm,n,v90) \n      call vector(n,v90,vnorm)\n      call orthog(vnorm,n,v90,v50)\n      call vector(n,v90,vnorm)\n      call normal(vnorm,n,v90) \n      call vector(n,v90,vnorm)\n      call orthog(vnorm,n,v90,v70)\n      call vector(n,v90,vnorm)\n      call normal(vnorm,n,v90) \n      call vector(n,v90,vnorm)\n\nc     orthog v90 to c60\n      call orthog(vnorm,n,c60,v90)\n      call vector(n,c60,vnorm)\n \n      call checko(n,v90,so)\n      call checko(n,v90,v10)\n      call checko(n,v90,v20)\n      call checko(n,v90,v30)\n      call checko(n,v90,v40)\n      call checko(n,v90,v50)\n      call checko(n,v90,v60)\n      call checko(n,v90,v70)\n      call checko(n,v90,v80)\n      call dotprod(k9,n,v90,c60)\n      write(9,*)'dotp v90 c60',k9\nc      call cv(v90,n)\n      write(12,*)'c60 dipole field after v90'\n      call output(c60,n)\n      call cv(c60,n)\n      call sumit(n,c60,sumo)\nc     ************************************************************\nc     v100\n      write(12,*)'v100'\n      do i=1,n\n          v100k=(230945*z(i)**9)-(437580*z(i)**7)+(270270*z(i)**5)\n          v100k=v100k-(60060*(z(i)**3))+(3465*(z(i)))\n          v100(i,1)=-(1.0d0\/128.0d0)*x(i)*z(i)*v100k\n          v100(i,2)=-(1.0d0\/128.0d0)*y(i)*z(i)*v100k\n          v100(i,3)=(1.0d0\/128.0d0)*(1.0d0-z(i)**2)*v100k\n      enddo\n      call output(v100,n)    \n      write(10,*)'v100 raw'\n      call checko(n,v100,so)\n      call checko(n,v100,v10)\n      call checko(n,v100,v20)\n      call checko(n,v100,v30)\n      call checko(n,v100,v40)\n      call checko(n,v100,v50)\n      call checko(n,v100,v60)\n      call checko(n,v100,v70)\n      call checko(n,v100,v80)\n      call checko(n,v100,v90)\n      write(10,*)'end raw'\n      \n      call dotprod(dotp,n,v100,c60)\n      write(9,*)'dotp v100 c60',dotp\nc     ************************************************************\nc     v110\n      write(12,*)'v110'\n      do i=1,n\n          v110k=(969969.0d0*(z(i)**10.0d0))-(2078505.0d0*(z(i)**8))\n          v110k=v110k+(1531530.0d0*(z(i)**6))-(450450.0d0*(z(i)**4))\n          v110k=v110k+(45045.0d0*(z(i)**2.0d0))-693.0d0\n          v110(i,1)=(-1.0d0\/256.0d0)*x(i)*z(i)*v110k\n          v110(i,2)=(-1.0d0\/256.0d0)*y(i)*z(i)*v110k\n          v110(i,3)=(1.0d0\/256.0d0)*(1.0d0-z(i)**2)*v110k\n      enddo\n      call output(v110,n)\n      write(12,*)'c60 upto 110'\n      call output(c60,n)\n      write(10,*)'v110 raw'\n      call checko(n,v110,so)\n      call checko(n,v110,v10)\n      call checko(n,v110,v20)\n      call checko(n,v110,v30)\n      call checko(n,v110,v40)\n      call checko(n,v110,v50)\n      call checko(n,v110,v60)\n      call checko(n,v110,v70)\n      call checko(n,v110,v80)\n      call checko(n,v110,v90)\n      call checko(n,v110,v100)\n      write(10,*)'end raw'\n \n      call orthog(vnorm,n,v110,so)\n      call vector(n,v110,vnorm)\n      call normal(vnorm,n,v110) \n      call vector(n,v110,vnorm)\n      call orthog(vnorm,n,v110,v10)\n      call vector(n,v110,vnorm)\n      call normal(vnorm,n,v110) \n      call vector(n,v110,vnorm)\n      call orthog(vnorm,n,v110,v30)\n      call vector(n,v110,vnorm)\n      call normal(vnorm,n,v110) \n      call vector(n,v110,vnorm)\n      call orthog(vnorm,n,v110,v50)\n      call vector(n,v110,vnorm)\n      call normal(vnorm,n,v110) \n      call vector(n,v110,vnorm)\n      call orthog(vnorm,n,v110,v70)\n      call vector(n,v110,vnorm)\n      call normal(vnorm,n,v110) \n      call vector(n,v110,vnorm)\n      call orthog(vnorm,n,v110,v90)\n      call vector(n,v110,vnorm)\n      call normal(vnorm,n,v110) \n      call vector(n,v110,vnorm)\n\nc     make orthog here\n      call orthog(vnorm,n,c60,v110)\n      call vector(n,c60,vnorm)\n\n      call checko(n,v110,so)\n      call checko(n,v110,v10)\n      call checko(n,v110,v20)\n      call checko(n,v110,v30)\n      call checko(n,v110,v40)\n      call checko(n,v110,v50)\n      call checko(n,v110,v60)\n      call checko(n,v110,v70)\n      call checko(n,v110,v80)\n      call checko(n,v110,v90)\n      call checko(n,v110,v100)\n      call dotprod(k11,n,v110,c60)\n      write(9,*)'dotp v110 c60',k11\nc      call cv(c60,n)\n      write(12,*)'c60 dipole field after v90'\n      call output(c60,n)\n      call cv(c60,n)\n      call sumit(n,c60,sumo)\nc     ***********************************************************\nc     v130\n      write(12,*)'v130'\n      do i=1,n\n          v130k=(16900975*(z(i)**12))-(44618574*(z(i)**10))\n          v130k=v130k+(43648605*(z(i)**8))-(19399380*(z(i)**6))\n          v130k=v130k+(3828825*(z(i)**4))-(270270*(z(i)**2))\n          v130k=v130k+3003\n          v130(i,1)=-(1.0d0\/1024.0d0)*x(i)*z(i)*v130k\n          v130(i,2)=-(1.0d0\/1024.0d0)*y(i)*z(i)*v130k\n          v130(i,3)=(1.0d0\/1024.0d0)*(1.0d0-z(i)**2)*v130k\n      enddo\n      call output(v130,n)\n      write(10,*)'v130 raw'\n      call checko(n,v130,so)\n      call checko(n,v130,v10)\n      call checko(n,v130,v20)\n      call checko(n,v130,v30)\n      call checko(n,v130,v40)\n      call checko(n,v130,v50)\n      call checko(n,v130,v60)\n      call checko(n,v130,v70)\n      call checko(n,v130,v80)\n      call checko(n,v130,v90)\n      call checko(n,v130,v100)\n      call checko(n,v130,v110)\n      write(10,*)'end raw'\n\nc     make orthog here\nc      call orthog(vnorm,n,c60,v130)\nc      call vector(n,c60,vnorm)\n  \n      call orthog(vnorm,n,v130,so)\n      call vector(n,v130,vnorm)\n      call normal(vnorm,n,v130) \n      call vector(n,v130,vnorm)\n      call orthog(vnorm,n,v130,v10)\n      call vector(n,v130,vnorm)\n      call normal(vnorm,n,v130) \n      call vector(n,v130,vnorm)\n      call orthog(vnorm,n,v130,v30)\n      call vector(n,v130,vnorm)\n      call normal(vnorm,n,v130) \n      call vector(n,v130,vnorm)\n      call orthog(vnorm,n,v130,v50)\n      call vector(n,v130,vnorm)\n      call normal(vnorm,n,v130) \n      call vector(n,v130,vnorm)\n      call orthog(vnorm,n,v130,v70)\n      call vector(n,v130,vnorm)\n      call normal(vnorm,n,v130) \n      call vector(n,v130,vnorm)\n      call orthog(vnorm,n,v130,v90)\n      call vector(n,v130,vnorm)\n      call normal(vnorm,n,v130) \n      call vector(n,v130,vnorm)\n      call orthog(vnorm,n,v130,v110)\n      call vector(n,v130,vnorm)\n      call normal(vnorm,n,v130) \n      call vector(n,v130,vnorm)\n\n \n      call checko(n,v130,so)\n      call checko(n,v130,v10)\n      call checko(n,v130,v20)\n      call checko(n,v130,v30)\n      call checko(n,v130,v40)\n      call checko(n,v130,v50)\n      call checko(n,v130,v60)\n      call checko(n,v130,v70)\n      call checko(n,v130,v80)\n      call checko(n,v130,v90)\n      call checko(n,v130,v100)\n      call checko(n,v130,v110)\n      call dotprod(k13,n,v130,c60)\n      write(9,*)'dotp v130 c60',k13\n      call cv(v130,n)\nc     ************************************************************\nc     v150 here\n      write(12,*)'v150'\n      do i=1,n\n          v150k=(145422675*(z(i)**14))-(456326325*(z(i)**12))\n          v150k=v150k+(557732175*(z(i)**10))-(334639305*(z(i)**8))\n          v150k=v150k+(101846745*(z(i)**6))-(14549535*(z(i)**4))\n          v150k=v150k+(765765*(z(i)**2))-(6435)\n          v150(i,1)=-(1.0d0\/2048.0d0)*x(i)*z(i)*v150k\n          v150(i,2)=-(1.0d0\/2048.0d0)*y(i)*z(i)*v150k\n          v150(i,3)=(1.0d0\/2048.0d0)*(1.0d0-z(i)**2)*v150k\n      enddo\n      call output(v150,n)\n      call normal(vnorm,n,v150)\n      call vector(n,v150,vnorm)\n      write(10,*)'v150 raw'\n      call checko(n,v150,so)\n      call checko(n,v150,v10)\nc     call checko(n,v150,v20)\n      call checko(n,v150,v30)\nc      call checko(n,v150,v40)\n      call checko(n,v150,v50)\nc      call checko(n,v150,v60)\n      call checko(n,v150,v70)\nc      call checko(n,v150,v80)\n      call checko(n,v150,v90)\nc      call checko(n,v150,v100)\n      call checko(n,v150,v110)\n      call checko(n,v150,v130)\n      write(10,*)'end raw'\n      call normal(vnorm,n,v150) \n      call vector(n,v150,vnorm)\n      \n      call orthog(vnorm,n,v150,so)\n      call vector(n,v150,vnorm)\n      call normal(vnorm,n,v150) \n      call vector(n,v150,vnorm)\n      call orthog(vnorm,n,v150,v10)\n      call vector(n,v150,vnorm)\n      call normal(vnorm,n,v150) \n      call vector(n,v150,vnorm)\n      call orthog(vnorm,n,v150,v20)\n      call vector(n,v150,vnorm)\n      call normal(vnorm,n,v150) \n      call vector(n,v150,vnorm)\n      call orthog(vnorm,n,v150,v30)\n      call vector(n,v150,vnorm)\n      call normal(vnorm,n,v150) \n      call vector(n,v150,vnorm)\n      call orthog(vnorm,n,v150,v50)\n      call vector(n,v150,vnorm)\n      call normal(vnorm,n,v150) \n      call vector(n,v150,vnorm)\n      call orthog(vnorm,n,v150,v70)\n      call vector(n,v150,vnorm)\n      call normal(vnorm,n,v150) \n      call vector(n,v150,vnorm)\n      call orthog(vnorm,n,v150,v90)\n      call vector(n,v150,vnorm)\n      call normal(vnorm,n,v150) \n      call vector(n,v150,vnorm)\n      call orthog(vnorm,n,v150,v110)\n      call vector(n,v150,vnorm)\n      call normal(vnorm,n,v150) \n      call vector(n,v150,vnorm)\n      call orthog(vnorm,n,v150,v130)\n      call vector(n,v150,vnorm)\n      call normal(vnorm,n,v150) \n      call vector(n,v150,vnorm)\n\n\n      call checko(n,v150,so)\n      call checko(n,v150,v10)\nc      call checko(n,v150,v20)\n      call checko(n,v150,v30)\nc      call checko(n,v150,v40)\n      call checko(n,v150,v50)\nc      call checko(n,v150,v60)\n      call checko(n,v150,v70)\nc      call checko(n,v150,v80)\n      call checko(n,v150,v90)\nc      call checko(n,v150,v100)\n      call checko(n,v150,v110)\n      call checko(n,v150,v130)\n \n      call dotprod(k15,n,v150,c60)\n      write(9,*)'dotp v150 c60',k15\n      call cv(v150,n)\nc     ************************************************************\nc     Calculations\nc     a(cage)=coeff(scalar).vsh sum off etc \n\n      do i=1,n\n         do j=1,3\n            cage(i,j)=(k0*so(i,j))+(k1*v10(i,j))\n            cage(i,j)=cage(i,j)+(k5*v50(i,j))+(k7*v70(i,j))\n            cage(i,j)=cage(i,j)+(k9*v90(i,j))\nc            cage(i,j)=cage(i,j)+(k11*v110(i,j))\nc            cage(i,j)=cage(i,j)+(k13*v130(i,j))\n            cage(i,j)=cage(i,j)+(k15*v150(i,j))\nc            cage(i,j)=cage(i,j)+(k17*v170(i,j))\n         enddo\n      enddo\nc   25 format(4x,d20.12,1x,d20.12,1x,d20.12)\n      do i=1,n\n         x(i)=cage(i,1)\n         y(i)=cage(i,2)\n         z(i)=cage(i,3)\n         write(14,25)x(i),y(i),z(i)\n      enddo\n\n      call orthog(vnorm,n,c60,cage)\n      call vector(n,c60,vnorm)\nc      call normal(vnorm,n,c60) \nc      call vector(n,c60,vnorm)\n      write(12,*)'c60 orthog cage'\n      call output(c60,n)\n\nc     testing how far away as same vector \n      call dotprod(dotp,n,c60,c60)\n      write(14,*)'dot c60 c60',dotp\n      call dotprod(dotp,n,c60,cage)\n      write(14,*)'dot c60 cage',dotp\nc     ************************************************************\n      stop\n      end\n      subroutine dotprod(dotp, n, a, b)\n      implicit real*8(a-h,o-z)\n      integer n\n      double precision dotp\n      dimension a(n,3), b(n,3)\n      dotp=0.0d0\n      do i=1,n\n         do j=1,3\n            dotp=dotp+(a(i,j)*b(i,j))\n         enddo\n      enddo\n      return\n      end\nc     __________________________________________________________\n      subroutine normal(vnorm,n,a)\n      implicit real*8(a-h,o-z)\n      integer n\n      double precision dotp,sqrk\n      dimension vnorm(n,3),a(n,3)\n      call dotprod(dotp,n,a,a)\n      do i=1,n\n         do j=1,3\n           sqrk=SQRT(abs(dotp))\n           vnorm(i,j)=a(i,j)\/sqrk\n         enddo\n      enddo\n      return\n      end\nc     __________________________________________________________\n      subroutine orthog(c,n,a,b)\n      implicit real*8(a-h,o-z)\n      integer n\n      double precision dotp\n      dimension vorth(n,3),a(n,3),b(n,3),c(n,3)\n      call dotprod(dotp,n,a,b)\n      do i=1,n\n         do j=1,3\n            vorth(i,j)=a(i,j)-(dotp*b(i,j))      \n         enddo\n      enddo\nc      call normal(c,n,vorth)\n      call vector(n,c,vorth)\n      return\n      end\nc     __________________________________________________________\n      subroutine output(vout, n)\n      implicit real*8(a-h,o-z)\n      integer n\n      dimension vout(n,3)\n   99 format(I2,1x,f24.19,1x,f24.19,1x,f24.19)\n      do i=1,n\n         write(12,99)i,vout(i,1), vout(i,2), vout(i,3)\n      enddo\n      return\n      end\nc     __________________________________________________________\n      subroutine checkn(vnorm, n)\n      implicit real*8(a-h,o-z)\n      integer n\n      double precision dotp\n      dimension vnorm(n,3)\n      call dotprod(dotp, n, vnorm, vnorm)\n      write(10,*)'Normalisation this should be  = 1',dotp\n      return\n      end\nc     _________________________________________________________\n      subroutine checko(n,vx,vy)\n      implicit real*8(a-h,o-z)\n      integer n\n      double precision dotp\n      dimension vx(n,3),vy(n,3)\n      call dotprod(dotp,n,vx,vy)\n      write(10,*)'Orthogalisation this should be  = 0',dotp\n      return\n      end\nc     _________________________________________________________\n      subroutine vector(n,a,b)\n      implicit real*8(a-h,o-z)\n      integer n\n      dimension a(n,3),b(n,3)\n      do i=1,n\n         do j=1,3\n           a(i,j)=b(i,j)\n         enddo\n      enddo\n      return\n      end\nc     _________________________________________________________\n      subroutine cv(x,n)\n      implicit real*8(a-h,o-z)\n      integer n\n      dimension x(n,3),sum(3)\n      do i=1,n\n         sum(1)=sum(1)+x(i,1)\n         sum(2)=sum(2)+x(i,2)\n         sum(3)=sum(3)+x(i,3)\n      enddo\n      write(12,*)'x',sum(1),' y',sum(2),' z',sum(3)\n      return\n      end\nc     ___________________________________________________\n      subroutine sumit(n,a,sumo)\n      implicit real*8(a-h,o-z)\n      integer n\n      double precision sumo\n      dimension a(n,3)\n      sumo=0.0d0\n      do i=1,n\n         do j=1,3\n         sumo=sumo+abs(a(i,j))\n         enddo\n      enddo\n      write(12,*)'sumo',sumo\n      return\n      end\n<\/pre>\n\n\n\n
<\/p>\n","protected":false},"excerpt":{"rendered":"
icorot.f c icosahedron rotator program 0.1aimplicit real*8(a-h,o-z)parameter(n=60)double precision m,m2,xtemp,ytemp,ztemp,theta,a,b,pdouble precision conv,switch,raddimension x(n),y(n),z(n)open(12,file=’c60.xyz’)open(13,file=’c60copy.xyz’)open(14,file=’c60ico.xyz’)switch=0.0d0p=20.90515744515d010 format(1x,f24.20,1x,f24.20,1x,f24.20)81 format(f24.20,1x,f24.20,1x,I10)do i=1,nread(12,10)x(i),y(i),z(i)enddorad=3.1415926540d0\/180.0d0do i=1,nx(i)=x(i)*rady(i)=y(i)*radz(i)=z(i)*radenddodo i=1,nwrite(13,10)x(i),y(i),z(i)enddodo i=1,9999999if (switch.eq.1) goto 91p=p*rad+1.0E-14*rada=(y(1)\/rad)*sin(p)+(z(1)\/rad)*cos(p)b=(y(2)\/rad)*sin(p)+(z(2)\/rad)*cos(p)conv=abs(a-b)p=p\/radc write(6,*)p,a,b,conv,iwrite(6,81)p,conv,iif (conv.lt.1.0E-15) switch=1enddo91 continuep=p*raddo i=1,nytemp=0.0d0ztemp=0.0d0ytemp=y(i)*cos(p)-z(i)*sin(p)ztemp=y(i)*sin(p)+z(i)*cos(p)y(i)=ytempz(i)=ztempenddowrite(13,*)’rotated’do i=1,nx(i)=x(i)\/rady(i)=y(i)\/radz(i)=z(i)\/radenddodo i=1,nwrite(13,10)x(i),y(i),z(i)write(14,10)x(i),y(i),z(i)enddostopend pent.f     c implicit real*8 (a-h,o-z) character*4 itex,itexh,itexb,ibas,ibash,istg integer p, s double precision aconst,r,tconst,zdf,test,pol,rr,dota,dotb double precision ratph, rstart, rstep, rpent, rhex,switch double precision conv parameter (n=5) […]<\/p>\n","protected":false},"author":1,"featured_media":0,"parent":0,"menu_order":0,"comment_status":"closed","ping_status":"closed","template":"","meta":{"footnotes":""},"class_list":["post-312","page","type-page","status-publish","hentry"],"_links":{"self":[{"href":"https:\/\/www.divergentfocus.co.uk\/index.php?rest_route=\/wp\/v2\/pages\/312","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.divergentfocus.co.uk\/index.php?rest_route=\/wp\/v2\/pages"}],"about":[{"href":"https:\/\/www.divergentfocus.co.uk\/index.php?rest_route=\/wp\/v2\/types\/page"}],"author":[{"embeddable":true,"href":"https:\/\/www.divergentfocus.co.uk\/index.php?rest_route=\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.divergentfocus.co.uk\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=312"}],"version-history":[{"count":0,"href":"https:\/\/www.divergentfocus.co.uk\/index.php?rest_route=\/wp\/v2\/pages\/312\/revisions"}],"wp:attachment":[{"href":"https:\/\/www.divergentfocus.co.uk\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=312"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}