*Deck OPolyGen
C
C Polygen: a program to generate spherical polyhedra with triangular
C faces by C.S. Pomelli (cris@ibm550.icqem.cnr.it)
C An equilater division algorithm is used.
C
      subroutine OPolyGen(Ptype,Pflag,TsAre,TsNum,XEN,YEN,ZEN,REN,
     $  TsEff,CV,JTR)
      implicit real*8 (d,g)
      Parameter (MxTs=5000)
      Parameter (MxVert = MxTs/2 +2)
      parameter (Icosa = 1)
      parameter (NVI=12, NTI=20, NEI=30)
      parameter (NVP=32, NTP=60, NEP=90)
      parameter (NVT= 4, NTT=4 , NET= 6)
      integer i,j,k,l,m,n,ii,jj,NV,NT0,NE0,NF
      integer NTPT,NVPT,Ptype,TsNum,TsEff,Pflag
      integer NFI,NFP,NFT,NDI,NDP,NDT
      integer trvo(60,3),treo(60,3),edo(90,2)
      integer JTR(MxTs,3)
      integer oldtr(100,100),ednew(90,100),trnew(60,100,100)
      integer tettrv(3,4), tettre(3,4), teted(2,6)
      integer icotrv(3,20),icotre(3,20),icoed(2,30)
      integer pentrv(3,60),pentre(3,60),pened(2,90)
      Real*8 tetve(3,4),icove(3,12),penve(3,32),dnorm,TsAre,pi,v1(3),
     $  v2(3),v3(3),alpha,beta,costheta,theta,cos1,cos2,sintheta,XEN,
     $  YEN,ZEN,REN,CV(MxVert,3)
C     data pi/3.141592653589793d0/
C Data for original polyhedra
C 1) Tetrahedron
      data teted/1,2,2,3,1,3,1,4,3,4,2,4/
      data tettre/1,2,3,1,6,4,3,5,4,2,5,6/
      data tettrv/1,2,3,1,2,4,1,3,4,2,3,4/
      data tetve/
     $  -.577350269D+00, .577350269D+00, .577350269D+00,
     $  .577350269D+00,-.577350269D+00, .577350269D+00,
     $  -.577350269D+00,-.577350269D+00,-.577350269D+00,
     $  .577350269D+00, .577350269D+00,-.577350269D+00/
C 2) Icosaedron
      data icoed/
     $  1, 2, 2, 3, 1, 3, 2, 6, 1, 6, 3, 4, 1, 4, 4, 5, 1, 5, 5, 6, 3,
     $  7, 2, 7, 6,11, 2,11, 7,11, 4, 8, 3, 8, 7, 8, 5, 9, 4, 9, 8, 9,
     $  6,10, 5,10, 9,10,10,11, 8,12, 7,12,11,12, 9,12,10,12/
      data icotre/
     $  1, 2, 3, 1, 4, 5, 3, 6, 7, 7, 8, 9, 9,10, 5, 2,11,12,
     $  4,13,14,12,15,14, 6,16,17,11,18,17, 8,19,20,16,21,20,10,22,23,
     $  19,24,23,22,25,13,18,26,27,15,28,27,21,29,26,24,30,29,25,28,30/
      data icotrv/
     $  1 , 2, 3, 1, 2, 6, 1, 3, 4, 1, 4, 5, 1, 5, 6, 2, 3, 7,
     $  2, 6,11, 2, 7,11, 3, 4, 8, 3, 7, 8, 4, 5, 9, 4, 8, 9, 5, 6,10,
     $  5, 9,10, 6,10,11, 7, 8,12, 7,11,12, 8, 9,12, 9,10,12,10,11,12/
      data icove/
     $  .000000000D+00,.000000000D+00, .100000000D+01,.276393202D+00,
     $  .850650808D+00, .447213595D+00,-.723606798D+00, .525731112D+00,
     $  .447213595D+00,-.723606798D+00,-.525731112D+00, .447213595D+00,
     $ .276393202D+00,-.850650808D+00, .447213595D+00,.894427191D+00,-
     $ .219071479D-15, .447213595D+00,-.276393202D+00, .850650808D+00,-
     $ .447213595D+00,-.894427191D+00, .109535740D-15,-.447213595D+00,-
     $ .276393202D+00,-.850650808D+00,-.447213595D+00,.723606798D+00,-
     $ .525731112D+00,-.447213595D+00,.723606798D+00, .525731112D+00,-
     $ .447213595D+00,.000000000D+00, .000000000D+00,-.100000000D+01/
C 3) Pentakisdodecahedron
      data pened/
     $  1 , 2, 2, 3, 1, 3, 2, 6, 1, 6, 3, 4, 1, 4, 4, 5, 1, 5, 5, 6, 3,
     $  8,2, 8, 6, 7, 2, 7, 7,12, 2,12, 8,12, 4, 9, 3, 9, 8,13, 3,13, 9,
     $  13, 5,10, 4,10, 9,14, 4,14,10,14, 6,11, 5,11,10,15, 5,15,11,15,
     $  7,16, 6,16,11,16,12,17, 7,17,16,17,12,18, 8,18,13,18,13,19, 9,
     $  19,14,19,14,20,10,20,15,20,15,21,11,21,16,21,17,22,12,22,18,22,
     $  18,23,13,23,19,23,19,24,14,24,20,24,20,25,15,25,21,25,17,26,16,
     $  26,21,26,22,27,17,27,26,27,22,28,18,28,23,28,23,29,19,29,24,29,
     $  24,30,20,30,25,30,25,31,21,31,26,31,27,28,28,29,29,30,30,31,27,
     $  31,28,32,27,32,31,32,29,32,30,32/
      data pentre/
     $  1, 2, 3, 1, 4, 5, 3, 6, 7, 7, 8, 9, 9,10, 5, 2,11,12, 4,13,14,
     $  14,15,16,12,17,16, 6,18,19,11,20,21,19,22,21, 8,23,24,18,25,26,
     $  24,27,26,10,28,29,23,30,31,29,32,31,13,33,34,28,35,34,15,36,37,
     $  33,38,37,17,39,40,20,41,40,22,42,43,25,44,43,27,45,46,30,47,46,
     $  32,48,49,35,50,49,36,51,52,39,53,52,41,54,55,42,56,55,44,57,58,
     $  45,59,58,47,60,61,48,62,61,38,63,64,50,65,64,51,66,67,63,68,67,
     $  53,69,70,54,71,70,56,72,73,57,74,73,59,75,76,60,77,76,62,78,79,
     $  65,80,79,66,81,69,71,82,72,74,83,75,77,84,78,68,85,80,81,86,87,
     $  85,88,87,82,89,86,83,90,89,84,88,90/
      data pentrv/
     $  1, 2, 3, 1, 2, 6, 1, 3, 4, 1, 4, 5, 1, 5, 6, 2, 3, 8, 2, 6, 7,
     $  2, 7,12, 2, 8,12, 3, 4, 9, 3, 8,13, 3, 9,13, 4, 5,10, 4, 9,14,
     $  4,10,14, 5, 6,11, 5,10,15, 5,11,15, 6, 7,16, 6,11,16, 7,12,17,
     $  7,16,17, 8,12,18, 8,13,18, 9,13,19, 9,14,19,10,14,20,10,15,20,
     $  11,15,21,11,16,21,12,17,22,12,18,22,13,18,23,13,19,23,14,19,24,
     $  14,20,24,15,20,25,15,21,25,16,17,26,16,21,26,17,22,27,17,26,27,
     $  18,22,28,18,23,28,19,23,29,19,24,29,20,24,30,20,25,30,21,25,31,
     $  21,26,31,22,27,28,23,28,29,24,29,30,25,30,31,26,27,31,27,28,32,
     $  27,31,32,28,29,32,29,30,32,30,31,32/
      data penve/
     $  .000000000D+00, .000000000D+00, .100000000D+01, .491123473D+00,
     $  .356822090D+00, .794654472D+00,-.187592474D+00, .577350269D+00,
     $  .794654472D+00,-.607061998D+00, .540159669D-09, .794654472D+00,-
     $  .187592475D+00,-.577350269D+00, .794654472D+00, .491123473D+00,-
     $  .356822091D+00, .794654472D+00, .894427191D+00, .000000000D+00,
     $  .447213595D+00, .276393203D+00, .850650808D+00, .447213595D+00,-
     $  .723606797D+00, .525731113D+00, .447213595D+00,-.723606799D+00,-
     $  .525731111D+00, .447213595D+00, .276393201D+00,-.850650809D+00,
     $  .447213595D+00, .794654472D+00, .577350269D+00, .187592474D+00,-
     $  .303530999D+00, .934172359D+00, .187592474D+00,-.982246946D+00,
     $  .873996703D-09, .187592474D+00,-.303531000D+00,-.934172359D+00,
     $  .187592474D+00, .794654471D+00,-.577350271D+00, .187592474D+00,
     $  .982246946D+00, .000000000D+00,-.187592474D+00, .303530999D+00,
     $  .934172359D+00,-.187592474D+00,-.794654472D+00, .577350270D+00,-
     $  .187592474D+00,-.794654473D+00,-.577350268D+00,-.187592474D+00,
     $  .303530997D+00,-.934172359D+00,-.187592474D+00, .723606798D+00,
     $  .525731112D+00,-.447213596D+00,-.276393202D+00, .850650808D+00,-
     $  .447213596D+00,-.894427191D+00, .795855278D-09,-.447213596D+00,-
     $  .276393203D+00,-.850650808D+00,-.447213596D+00, .723606797D+00,-
     $  .525731113D+00,-.447213596D+00, .607061998D+00, .000000000D+00,-
     $  .794654472D+00, .187592474D+00, .577350269D+00,-.794654472D+00,-
     $  .491123473D+00, .356822090D+00,-.794654472D+00,-.491123474D+00,-
     $  .356822089D+00,-.794654472D+00, .187592473D+00,-.577350269D+00,-
     $  .794654472D+00,.000000000D+00, .000000000D+00,-.100000000D+01/
C
C set pi
C
      PI=3.14159265358979323846264338327950288419716939937510D+00
C
C  A seconda delle due opzioni di funzionamento (Area ottimale o
C  numero di tessere ottimale ) vengono stabilite il tipo di poliedro
C  (Ptype) e la  frequenza di divisione (NF)
C
C  Pflag = 0 Numero di tessere: si genera il poliedro disponibile
C  con il numero di tessere piu' prossimo a TsNum
C
C  Pflag  =1  Area: si genera il poliedro con il numero di tessere tali
C  da avere un area il piu' simile possibile
      if (Pflag.eq.1) then
        TsNum = INT (4.0d0 * pi * REN**2 / TsAre + 0.5 )
      endif
C
C  -si controlla che TsNum <= MxTs
C
      if (TsNum.gt.MxTs) Call GauErr('Too many tesserae in OPolyGen')
C
C  -calcola il valore di NF ottimale per ciascuna famiglia di poliedri.
C
      NFI = IGFix( SQRT(TsNum / ( NTI * 1.0d0 )) + 0.5d0)
      NFP = IGFix( SQRT(TsNum / ( NTP * 1.0d0 )) + 0.5d0)
      NFT = IGFix( SQRT(TsNum / ( NTT * 1.0d0 )) + 0.5d0)
C
C  -trova lo scostamento del numero di tessere in input
C
      NDI = IABS( TsNum - NTI * NFI**2 )
      NDP = IABS( TsNum - NTP * NFP**2 )
      NDT = IABS( TsNum - NTT * NFT**2 )
C
C  -sceglie tipo e frequenza di divisione
C  a parita' di differenza in numero di tessere sceglie quello con NF
C  minore
C
      if     ((NDP.le.NDI).and.(NDP.le.NDT)) then
        Ptype = 2
        NF = NFP
        TsEff = NTP * NFP**2
      else
        if (NDI.le.NDT) then
          Ptype = Icosa
          NF = NFI
          TsEff = NTI * NFI**2
        else
          Ptype = 3
          NF = NFT
          TsEff = NTT * NFT**2
        endif
      endif
C
C  -carica gli opportuni valori iniziali
C
      if (Ptype.eq.1) then
C
C  icosaedro
C
        NT0=NTI
        NV=NVI
        NE0=NEI
        do 2000 i=1, NT0
          do 2010 k=1,3
            trvo(i,k)=icotrv(k,i)
            treo(i,k)=icotre(k,i)
 2010 continue
 2000 continue
        do 2020 i=1,NV
          do 2030 k=1,3
            CV(i,k)=icove(k,i)
 2030 continue
 2020 continue
        do 2040 i=1,NE0
          do 2050 k=1,2
            edo(i,k)=icoed(k,i)
 2050 continue
 2040 continue
      elseif (Ptype.eq.2) then
C
C  -pentakisdodecaedro
C
        NT0=NTP
        NV=NVP
        NE0=NEP
        do 2060 i=1, NT0
          do 2070 k=1,3
            trvo(i,k)=pentrv(k,i)
            treo(i,k)=pentre(k,i)
 2070 continue
 2060 continue
        do 2080 i=1,NV
          do 2090 k=1,3
            CV(i,k)=penve(k,i)
 2090 continue
 2080 continue
        do 2100 i=1,NE0
          do 2110 k=1,2
            edo(i,k)=pened(k,i)
 2110 continue
 2100 continue
      elseif (Ptype.eq.3) then
C
C  tetraedro
C
        NT0=NTT
        NV=NVT
        NE0=NET
        do 2120 i=1, NT0
          do 2130 k=1,3
            trvo(i,k)=tettrv(k,i)
            treo(i,k)=tettre(k,i)
 2130 continue
 2120 continue
        do 2140 i=1,NV
          do 2150 k=1,3
            CV(i,k)=tetve(k,i)
 2150 continue
 2140 continue
        do 2160 i=1,NE0
          do 2170 k=1,2
            edo(i,k)=teted(k,i)
 2170 continue
 2160 continue
      endif
      NVPT=NV+1
C
C  -nuovi vertici posti lungo i vecchi spigoli
C  -le regole di calcolo derivano da  calcoli di algebra
C   lineare e trigonometria
C  -i vertici vengono memorizzati in ordine progressivo ed riferiti
C  tramite ednew al vertice di appartenenza
C
      do 2180 j=1,NE0
        do 2190 k=1,3
          v1(k)=CV(edo(j,1),k)
          v2(k)=CV(edo(j,2),k)
 2190 continue
        costheta=
     $    ( v1(1)*v2(1) + v1(2)*v2(2) + v1(3)*v2(3)) /
     $    (sqrt(v1(1)**2 + v1(2)**2 + v1(3)**2) *
     $     sqrt(v2(1)**2 + v2(2)**2 + v2(3)**2))
        theta=acos(costheta)
        sintheta=sin(theta)
        do 2200 l=1,NF-1
          m=NF-l
          cos1=cos(theta*l/NF)
          cos2=cos(theta*(NF-l)/NF)
          alpha=(cos1-costheta*cos2)/sintheta**2
          beta=(cos2-costheta*cos1)/sintheta**2
          do 2210 k=1,3
            v3(k)=alpha*v1(k)+beta*v2(k)
 2210 continue
          dnorm = sqrt(v3(1)**2 + v3(2)**2 + v3(3)**2)
          do 2220 k=1,3
            v3(k) = v3(k) /dnorm
 2220 continue
          do 2230 k=1,3
            CV(NVPT,k)=v3(k)
 2230 continue
          ednew(j,l+1)=NVPT
          NVPT=NVPT+1
 2200 continue
 2180 continue
C
C  -nuovi vertici non posti lungo i vecchi spigoli
C  -a partire dai vertici in ednew secondo regole analoghe alle
C  precedenti, vengono memorizzati a seconda del triangolo in trnew
C  trnew(triangolo,fila,n ordine)=nvertice
C
      do 2240 j=1,NT0
        ii=treo(j,1)
        jj=treo(j,3)
        do 2250 l=3,NF
          do 2260 k=1,3
            v1(k)=CV(ednew(ii,l),k)
            v2(k)=CV(ednew(jj,l),k)
 2260 continue
        costheta=
     $    ( v1(1)*v2(1) + v1(2)*v2(2) + v1(3)*v2(3)) /
     $    (sqrt(v1(1)**2 + v1(2)**2 + v1(3)**2) *
     $     sqrt(v2(1)**2 + v2(2)**2 + v2(3)**2))
          theta=acos(costheta)
          sintheta=sin(theta)
          do 2270 n=1,l-2
            cos1=cos(theta*n/(l-1))
            cos2=cos(theta*(l-1-n)/(l-1))
            alpha=(cos1-costheta*cos2)/sintheta**2
            beta=(cos2-costheta*cos1)/sintheta**2
            do 2280 k=1,3
              v3(k)=alpha*v1(k)+beta*v2(k)
 2280 continue
            dnorm = sqrt(v3(1)**2 + v3(2)**2 + v3(3)**2)
            do 2290 k=1,3
              v3(k) = v3(k) /dnorm
 2290 continue
            do 2300 k=1,3
              CV(NVPT,k)=v3(k)
 2300 continue
            trnew(j,l,n+1)=NVPT
            NVPT=NVPT+1
 2270 continue
 2250 continue
 2240 continue
      NV=NVPT-1
C
C  -ora per ogni triangolo originario vengono posti nella matrice oldtr
C  i numeri di ordine dei vertici originali,creati lungo gli spigoli e
C  no dei vecchi triangoli secondo lo schema:
C
C          11                    La disposizione e' secondo
C          | \                   la posizione geometrica del
C          |  \                  vertice.
C          21--22
C          | \  |\               i vertici (1,1)-(NF+1,1)-(NF+1,NF+1) sono
C          |  \ | \              quelli originari.
C          31--32--33
C          | \  |  | \           i nuovi triangoli sono:
C          |  \ |  |  \          (i,j)-(i+1,j)-(i+1,j+1)
C          41--42--43--44        i=1,...,NF j=1,...,i
C
C                                (i,j)--(i,j+1),(i+1,j+1)
C                                i=2,...,NF j=1,...,i-1
C
      NTPT=1
      do 2310 n=1,NT0
C
C  -1 vecchi spigoli
C
        oldtr(1,1)=trvo(n,1)
        oldtr(NF+1,1)=trvo(n,2)
        oldtr(NF+1,NF+1)=trvo(n,3)
C
C  -2 nuovi vertici lungo i vecchi spigoli
C
        do 2320 l=2,NF
          oldtr(l,1)=ednew(treo(n,1),l)
          oldtr(NF+1,l)=ednew(treo(n,2),l)
          oldtr(l,l)=ednew(treo(n,3),l)
 2320 continue
C
C  -3 nuovi vertici non lungo i vecchi spigoli
C
        do 2330 l=3,NF
          do 2340 m=2,l-1
            oldtr(l,m)=trnew(n,l,m)
 2340 continue
 2330 continue
C
C  -ora si creano i nuovi triangoli
C
        do 2350 i=1,NF
          do 2360 j=1,i
            JTR(NTPT,1)=oldtr(i,j)
            JTR(NTPT,2)=oldtr(i+1,j)
            JTR(NTPT,3)=oldtr(i+1,j+1)
            NTPT=NTPT+1
 2360 continue
 2350 continue
        do 2370 i=2,NF
          do 2380 j=1,i-1
            JTR(NTPT,1)=oldtr(i,j)
            JTR(NTPT,2)=oldtr(i,j+1)
            JTR(NTPT,3)=oldtr(i+1,j+1)
            NTPT=NTPT+1
 2380 continue
 2370 continue
 2310 continue
C
C scrittura della JVT
C
      do 2390 i=1,NV
         CV(i,1) = CV(i,1) * REN + XEN
         CV(i,2) = CV(i,2) * REN + YEN
         CV(i,3) = CV(i,3) * REN + ZEN
 2390 continue
      return
      end
