first commit

README.md

@@ -0,0 +1 @@
+# Geometry_and_Algorithms_for_CAD

cdg0gv-2010/beispiele/Makefile

@@ -0,0 +1,27 @@
+#PC=	ppc386
+PC=	ppcx64
+
+CDGDIR=         ..
+INCLDIR=        $(CDGDIR)/include
+UNITDIR=        $(CDGDIR)/units
+TOOLDIR=        $(CDGDIR)/tools
+
+PFLAGS=		-So -Fu../units -I../include
+
+%:		%.p
+		$(PC) $(PFLAGS)  $<
+
+units:		
+		(cd $(UNITDIR) && make units)
+		(cd $(TOOLDIR) && make tools)
+
+clean:
+		rm -f *~ *.o  
+
+cleanall:
+		rm -f *~ *.o n_eck 
+		(cd $(TOOLDIR) && make clean)
+		(cd $(UNITDIR) && make clean)
+
+
+

cdg0gv-2010/beispiele/README.fpc

@@ -0,0 +1,16 @@
+Installationschritte:
+
+1. "make" ausfuehren
+   erzeugt 
+   a) die notwendigen units in Dir UNITS und
+   b) pldv in TOOLS 
+
+2. Compilieren des Programms n_eck.p:  "make n_eck"
+   Ausfuehren von n_eck: "n_eck"
+
+3. "make clean" loescht .o Dateien und n_eck und tori_h
+4. "make cleanall" loescht 
+    a) ALLE .o-Dateien und 
+    b) .o- und .ppu-Dateien in UNITS
+    c) pldv in TOOLS 
+

cdg0gv-2010/beispiele/bezkur.p

@@ -0,0 +1,139 @@
+program bezkur;
+ uses geograph;
+ const n_bez=20;
+ type  bez_array   = array[0..n_bez] of real; 
+       bezpt_array = array[0..n_bez] of vt2d;
+ var   p : vts2d;         bx,by,bxnew,bynew : bez_array;
+       dt,t : real;       bezpts: bezpt_array;   curvept: vt2d;
+       i,j,degree,ipld,np,ivorn,degree1,nelev : integer;
+{******}
+
+function bezier_comp(degree: integer; var coeff: bez_array; t: real) : real;
+{Berechnet eine Komponente einer Bezier-Kurve.}
+{aus FARIN: Curves and surfaces...}
+ var i,n_choose_i : integer;    fact,t1,aux : real;
+ begin
+   t1:= 1-t;  fact:=1;  n_choose_i:= 1;
+   aux:= coeff[0]*t1;
+   for i:= 1 to degree-1 do
+       begin
+       fact:= fact*t;
+       n_choose_i:= n_choose_i*(degree-i+1) div i;
+       aux:= (aux + fact*n_choose_i*coeff[i])*t1;
+       end;
+   aux:= aux + fact*t*coeff[degree];
+   bezier_comp:= aux;
+ end;  {bezier_comp}
+{*************}
+
+function bezier_comp_decas(degree: integer; var coeff: bez_array; t: real) : real;
+{Berechnet eine Komponente einer Bezier-Kurve mit DECASTELJAU_Alg..}
+{aus FARIN: Curves and Surfaces, S. 33,34}
+ var i,r : integer;    t1 : real;       coeffa: bez_array;
+ begin
+   t1:= 1-t;  
+   coeffa:= coeff;
+   for r:= 1 to degree do
+       for i:= 0 to degree-r do coeffa[i]:= t1*coeffa[i] + t*coeffa[i+1];
+   bezier_comp_decas:= coeffa[0];
+ end;  {bezier_comp_decas}
+{*************}
+
+procedure degree_elev_comp(degree: integer; var coeff: bez_array; 
+                                            var coeffnew: bez_array);
+{Berechnet die Koeffizienten f"ur eine Graderh"ohung.}
+var i,degree1: integer; 
+begin
+  degree1:= degree+1;
+  coeffnew[0]:= coeff[0];  coeffnew[degree1]:= coeff[degree];
+  for i:= 1 to degree do   
+      coeffnew[i]:= (i*coeff[i-1] + (degree1-i)*coeff[i])/degree1;
+end;  { degree_elev_comp }
+{************}
+
+procedure beziercurvept2d_decas(degree: integer; var bezpts: bezpt_array; t: real; 
+                                                 var curvept: vt2d); 
+{Berechnet und zeichnet einen Punkt einer Bezier-Kurve mit DECASTELJAU_Alg..}
+{aus FARIN: Curves and Surfaces, S. 33,34}
+ var i,r : integer;    t1 : real;       bezptsa: bezpt_array;
+ begin
+   t1:= 1-t;  
+   bezptsa:= bezpts;
+   for r:= 1 to degree do
+       for i:= 0 to degree-r do 
+        begin
+        lcomb2vt2d(t1,bezptsa[i],t,bezptsa[i+1], bezptsa[i]);
+        point2d(bezptsa[i],0); 
+        if (i>0) and (r<degree) then line2d(bezptsa[i-1],bezptsa[i],0);
+        end;
+   curvept:= bezptsa[0];
+ end;  {beziercurvept2d_decas}
+{*************}
+
+{*************}
+ begin {Hauptprogramm}
+   writeln('pld-Datei ? (ja:1)');               read(ipld);
+   graph_on(ipld);
+   repeat
+     writeln('*** ebene BEZIER-Kurve (mit CASTELJAU oder hornbez) *** ');
+     writeln('Welcher Grad ? ');                  readln(degree);
+     writeln('Kontroll-Punkte (xi,zi):  ');
+     for i:= 0 to degree do
+         begin
+         write(i,'. Punkt ? ');                   readln(bx[i],by[i]);
+         end;
+     writeln('Anzahl der Kurven-Punkte:  np ?');  readln(np);
+     dt:= 1/np;     t:= 0;
+     for i:= 0 to np do
+         begin
+         p[i].x:= bezier_comp_decas(degree,bx,t);
+         p[i].y:= bezier_comp_decas(degree,by,t);
+         t:= t+dt;
+         end;
+     draw_area(180,180,20,20,1.3);
+     for i:= 1 to degree do
+         begin
+         pointc2d(bx[i],by[i],0);
+         linec2d(bx[i-1],by[i-1],bx[i],by[i],1);
+         end;
+     pointc2d(bx[0],by[0],0);
+     new_linewidth(2);
+     curve2d(p,0,np,0);   
+     new_linewidth(1);
+     {arrowc2d(0,0,60,0,2);  arrowc2d(0,0,0,60,2);}
+(*     
+{ Casteljau-Alg. fuer t=...:}
+     for i:= 0 to degree do put2d(bx[i],by[i], bezpts[i]);
+     t:= 0.6;
+     beziercurvept2d_decas(degree,bezpts,t,curvept); 
+*)
+
+{Graderhoehung:}
+     nelev:= 10;
+     for j:= 1 to nelev do
+        begin
+        degree_elev_comp(degree,bx,bxnew);
+        degree_elev_comp(degree,by,bynew);
+        degree1:= degree+1;
+        if j=nelev then
+        for i:= 0 to degree1 do pointc2d(bxnew[i],bynew[i],0);
+        if j=nelev then 
+        for i:= 0 to degree1-1 do linec2d(bxnew[i],bynew[i],bxnew[i+1],bynew[i+1],0);
+        bx:= bxnew;  by:= bynew;  degree:= degree1;
+        end;
+
+     draw_end;
+     writeln(' von vorn ? (nein: 0) ');          readln(ivorn);
+   until ivorn=0;
+   graph_off;
+ end.
+
+
+
+
+
+
+
+
+
+

cdg0gv-2010/beispiele/cassini.p

@@ -0,0 +1,150 @@
+program cassini;
+{mit Methode des "steilsten Weges"}
+uses linux,geograph;
+type  
+      funct2d = function(x,y: real) : real;    { fuer implicit_curvepts }
+      funct2d2d = procedure(x,y: real; var gx,gy: real);
+
+var   startpt,endpt,start_tangentvt: vt2d;
+      tangent_stepl,curvept_delta : real;
+      n2max,n2 : integer;
+      f: funct2d;
+      gradf: funct2d2d;
+      p : vts2d;
+
+ var   a,c,scalefac,x0,y0 : real;
+       p0,p1,p2 : vt2d;
+       inz,i,ipld : integer;
+{******}
+
+procedure implicit_curvepts(startpt,endpt,start_tangentvt: vt2d;
+                             tangent_stepl,curvept_delta : real;
+                                                   n2max : integer;
+                                                        f: funct2d;
+                                                    gradf: funct2d2d;
+                                                   var n2: integer;
+                                                   var p : vts2d);
+{berechnet Punkte p[0],...,p[n2] einer impliziten Kurve f(x,y)=0 mit der
+ Methode des "steilsten Weges".
+ startpt=Startpunkt, endpt=Endpunkt,
+ tangent_stepl= Schrittweite in Tangentenrichtung,
+ curvept_delta= Genauigkeitsschranke zur Berechnung von Kurvenpunkten,
+ n2max= maximales n2, n2max muss kleiner als array_size sein !.}
+ var i : integer;    tv,dv,ps,curve_endpt : vt2d;
+     g,ghn,gx,gy,fac,test,stepl,delta,cc : real;
+{****}
+procedure newton_step(xi,yi : real; var xi1,yi1 :real);
+ var  fi,gxi,gyi,t,cc : real;
+ begin
+   fi:= f(xi,yi);            gradf(xi,yi, gxi,gyi);
+   cc:= sqr(gxi) + sqr(gyi);
+   if cc>0 then t:= -fi/cc else begin t:= 0; pointc2d(xi,yi,0); end;
+   xi1:= xi + t*gxi;  yi1:= yi + t*gyi;
+ end;
+{***}
+procedure curve_point(p0 : vt2d; var pc : vt2d);
+{Sucht einen Punkt der Kurve entlang des steilsten Weges}
+ var xi,yi,xi1,yi1,delta : real;
+ begin
+   xi:= p0.x;  yi:= p0.y;
+   repeat                                      {pointc2d(xi,yi,10);}
+     newton_step(xi,yi,xi1,yi1);
+     delta:= abs(xi-xi1) + abs(yi-yi1);
+     xi:= xi1;  yi:= yi1;
+   until delta < curvept_delta;
+   put2d(xi1,yi1, pc);
+ end;   { curve_point }
+{*********}
+ begin {implicitcurve}      {Mit while-Schleife, n2 ist "var"}
+   curve_point(startpt, p[0]);       curve_point(endpt, curve_endpt);
+   stepl:= tangent_stepl;            delta:= stepl;
+   tv:= start_tangentvt;             fac:= stepl/length2d(tv);
+   lcomb2vt2d(1,p[0], fac,tv, ps);   curve_point(ps,p[1]);
+   i:=1;   test:= 0;
+   while (delta>0.7*stepl) and (i<n2max) do
+       begin
+       diff2d(p[i],p[i-1], dv);
+       gradf(p[i].x,p[i].y, gx,gy);   put2d(-gy,gx, tv);
+       cc:= length2d(tv);
+       if cc>0 then
+          begin
+          fac:= stepl/cc;
+          test:= scalarp2d(tv,dv);       if test<0 then fac:= -fac;
+          lcomb2vt2d(1,p[i], fac,tv, ps);
+          end
+          else lcomb2vt2d(2,p[i],1,p[i-1], ps);
+       i:= i+1;
+       curve_point(ps,p[i]);
+       delta:= distance2d(p[i],curve_endpt);
+       end;  { while }
+   if i<n2max then begin n2:= i+1; p[n2]:= curve_endpt; end
+              else n2:= i;
+  end; {implicit_curvepts}
+{***************}
+
+function fcass(x,y : real) : real;
+ begin
+   fcass:= sqr(x*x+y*y) - 2*c*c*(x*x-y*y) - (a*a*a*a-c*c*c*c) ;
+ end;
+{*****}
+procedure gradfcass(x,y : real; var gx,gy: real);
+ begin
+   gx:= 4*x*(x*x+y*y) - 4*c*c*x;
+   gy:= 4*y*(x*x+y*y) + 4*c*c*y;
+ end;
+{*****}
+
+{***********************}
+ begin {Hauptprogramm}
+   writeln('pld-Datei ? (ja:1, nein:0)');                 readln(ipld);
+   graph_on(ipld);
+   writeln('**************************************************');
+   writeln('            Implizite Kurven mit ');
+   writeln(        'Methode des steilsten Weges');
+   writeln('**************************************************');
+   repeat
+     writeln(' CASSINIsche Kurve (Bronstein S.141):');
+     writeln('(x*x+y*y)**2 -2*c*c*(x*x-y*y) - (a**4-c**4) = 0 ');
+     writeln('Fuer a=c ergibt sich eine LEMNISKATE');
+     writeln('  a ?  (>0) ,  c (0.5 ... 2) ?');          readln(a,c);
+{     writeln('Anfangspunkt (x0,y0) ? (nicht (0,0) !)');  readln(x0,y0);}
+     x0:= 1; y0:= 1; 
+     put2d(x0,y0, p0);     scalefac:= 30;
+{Zeichenflaeche :  }
+     draw_area(180,140,90,65,scalefac);
+     put2d(-2,0, p1);  put2d(2,0, p2);     arrow2d(p1,p2,2);
+     put2d(0,-1, p1);  put2d(0,1, p2);     arrow2d(p1,p2,2);
+
+     put2d(x0,y0,startpt); new_color(red); point2d(startpt,0); new_color(default);
+     endpt:= startpt;  put2d(0,1,start_tangentvt);
+     tangent_stepl:= 0.05;  curvept_delta:= eps6;  n2max:= 300;
+     f:= fcass;  gradf:= gradfcass;
+     implicit_curvepts(startpt,endpt,start_tangentvt,
+                             tangent_stepl,curvept_delta,
+                             n2max,f,gradf,n2,p);
+writeln('Anzahl der Punkte: ',n2);
+     curve2d(p,0,n2,0);
+     if c>a then
+        begin
+        for i:= 0 to n2 do p[i].x:= -p[i].x;  curve2d(p,0,n2,0);
+        end;
+     draw_end;          writeln;
+     writeln(' Noch eine Zeichnung ? (nein: 0) ');     read(inz);
+   until inz=0;
+   graph_off;
+ end.
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

cdg0gv-2010/beispiele/flaech_h.p

@@ -0,0 +1,153 @@
+{********************************************************************}
+{***  Projektion einer param. Flaeche mit UP cp_lines_before_convex_faces  ***}
+{********************************************************************}
+program flaech_h;
+uses geograph,hiddenl;
+(*  {const type und var in unit hiddenl:}
+const {Achtung: es muss  array_size>=nfmax  sein  !!!}
+       nfmax= 10000;  nemax=20000;    nsegmax=10;  npfmax=10;
+
+type   vts2d_pol  = array[0..npfmax] of vt2d;
+       vts3d_pol  = array[0..npfmax] of vt3d;
+       r_array_seg  = array[0..nsegmax] of real;
+       i_array_seg  = array[0..nsegmax] of integer;
+       box3d_dat = record
+                   xmin,xmax,ymin,ymax,zmin,zmax : real;
+                   end;
+       face_dat = record
+                  npf,nef : integer;
+                  fp,fe : array[1..npfmax] of integer;
+                  vis : boolean;
+                  box : box3d_dat;
+                  discentre,d : real;
+                  nv : vt3d;
+                  end;
+       edge_dat = record
+                  vis : boolean;
+                  ep1,ep2,color,linew : integer;
+                  end;
+
+var    ne,nf,np: integer;  {Anzahl der Kanten, Facetten,Punkte}
+       p    : vts3d;       {Punkte des Polyeders} 
+       face : array[1..nfmax] of face_dat;
+       edge : array[1..nemax] of edge_dat;
+       pdist: r_array;     {pdist[i]: Abstand d. i-ten Punktes von d. Bildeb.} 
+       error,oriented_faces,is_permitted,newstyles: boolean;
+*)
+
+ var   u,v,du,dv,u1,u2,v1,v2,r1 : real;
+       n1,n2,i,k,ik,iachs,inz,ianf,iplot,nf1,ne0 : integer;
+       ps1,ps2			: vt3d;     {fuer Schnittkanten der Flaechen}
+       inters			: boolean;
+{****************}
+
+begin {Hauptprogramm}
+   writeln('PLD-Datei ? (ja: 1)');  readln(iplot);
+   graph_on(iplot);
+   repeat
+     writeln('***  Projektion parametrisierter Flaechen         ***');    writeln;
+     writeln('***  hier: "Affensattel" und Zylinder             ***');
+      ne0:= 0;
+{--------------------------------}
+{ Affensattel: Typ ist "quadrangle"}
+     n1:= 30;    n2:= 30;   {Anzahl der Unterteilungen im Parameterbereich};
+     u1:= -0.9;  u2:= 0.9;  {Parametergrenzen}
+     v1:= -0.9;  v2:= 0.9;
+     du:= (u2-u1)/(n1-1);     dv:= (v2-v1)/(n2-1);
+     v:= v1;
+     for k:= 1 to n2 do
+         begin
+         u:= u1;
+         for i:= 1 to n1 do
+             begin
+             ik:=i + (k-1)*n1;
+             put3d(u,v,u*u*u-3*u*v*v, p[ik]);    {Parameterdarstellung}
+             u:= u+du;
+             end;
+         v:= v+dv;
+         end;
+     np:= n1*n2;
+{     aux_quadrangle_triang(n1,n2,true);}
+    aux_quadrangle(n1,n2,0,0,0); 
+     writeln(' np:',np);     writeln(' nf:',nf);     writeln(' ne:',ne);
+     nf1:= nf;
+     for i:= ne0+1 to ne do
+	  with edge[i] do begin color:=cyan; linew:=1; end;
+      ne0:= ne;
+{--------------------------------}      
+{Zylinder 1: x**2 + z**2 -r1**2=0, Typ ist "cylinder"}
+     r1:= 0.6;
+     n1:= 50;    n2:= 10;
+     v1:= -1.5;  v2:= 1;
+     u1:= 0;     u2:= pi2;
+     du:= (u2-u1)/n1;     dv:= (v2-v1)/(n2-1);
+     v:= v1;
+     for k:= 1 to n2 do
+         begin
+         u:= u1;
+         for i:= 1 to n1 do
+             begin
+             ik:=i + (k-1)*n1;
+             put3d(r1*cos(u),r1*sin(u),v, p[np+ik]);
+             u:= u+du;
+             end;
+         v:= v+dv;
+         end;
+     aux_cylinder(n1,n2,np,ne,nf);
+     for i:= ne0+1 to ne do
+	  with edge[i] do begin color:=blue; linew:=1; end;
+      ne0:= ne;
+{---------------------------------------------}
+
+{Durchdringung der beiden Flaechen:}
+    boxes_of_faces;
+    for i:= 1 to nf1 do         { 1.Flaechenschleife }
+        begin
+        for k:= nf1+1 to nf do  { 2.Flaechenschleife }
+            begin
+            is_face_face(i,k, ps1,ps2,inters);
+            if inters then
+               begin
+               p[np+1]:= ps1;  p[np+2]:= ps2;
+               with edge[ne+1] do 
+                    begin ep1:= np+1; ep2:= np+2; vis:= true; end;
+               np:= np+2;  ne:= ne+1;
+               with face[i] do begin fe[nef+1]:= ne; nef:= nef+1;  end;
+               with face[k] do begin fe[nef+1]:= ne; nef:= nef+1;  end;
+               end;  { if }
+            end;  { for k }
+        end;  { for i }
+     for i:= ne0+1 to ne do
+	  with edge[i] do begin color:=red; linew:=1; end;
+      ne0:= ne;
+{---------------------------------------------}
+      
+     repeat
+       writeln;
+       init_centralparallel_projection(2);
+       writeln(' Koordinaten-Achsen ? (Ja = 1)');
+       readln(iachs);
+{ Zeichnen : }
+       draw_area(200,200,100,100,50);
+       if iachs=1 then begin  cp_axes(1.8);  point2d(null2d,0);  end;
+
+       oriented_faces:= false;
+       is_permitted:= true;
+	newstyles:= true;
+       cp_lines_before_convex_faces(oriented_faces,is_permitted,newstyles);
+
+       draw_end ;        writeln ;
+       writeln(' Noch eine Projektion ? ( Ja = 1 )');
+     readln(inz);
+     until inz=0;
+
+   writeln('Noch einmal mit ANDEREN  n1, n2 ?  (ja: 1)');
+   readln(ianf);
+   until ianf=0;
+   graph_off;
+ end.
+
+
+
+
+

cdg0gv-2010/beispiele/flaech_h_off.p

@@ -0,0 +1,200 @@
+{********************************************************************}
+{***  Projektion einer param. Flaeche mit UP cp_lines_before_convex_faces  ***}
+{********************************************************************}
+program flaech_h;
+uses geograph,hiddenl;
+(*  {const type und var in unit hiddenl:}
+const {Achtung: es muss  array_size>=nfmax  sein  !!!}
+       nfmax= 10000;  nemax=20000;    nsegmax=10;  npfmax=10;
+
+type   vts2d_pol  = array[0..npfmax] of vt2d;
+       vts3d_pol  = array[0..npfmax] of vt3d;
+       r_array_seg  = array[0..nsegmax] of real;
+       i_array_seg  = array[0..nsegmax] of integer;
+       box3d_dat = record
+                   xmin,xmax,ymin,ymax,zmin,zmax : real;
+                   end;
+       face_dat = record
+                  npf,nef : integer;
+                  fp,fe : array[1..npfmax] of integer;
+                  vis : boolean;
+                  box : box3d_dat;
+                  discentre,d : real;
+                  nv : vt3d;
+                  end;
+       edge_dat = record
+                  vis : boolean;
+                  ep1,ep2,color,linew : integer;
+                  end;
+
+var    ne,nf,np: integer;  {Anzahl der Kanten, Facetten,Punkte}
+       p    : vts3d;       {Punkte des Polyeders} 
+       face : array[1..nfmax] of face_dat;
+       edge : array[1..nemax] of edge_dat;
+       pdist: r_array;     {pdist[i]: Abstand d. i-ten Punktes von d. Bildeb.} 
+       error,oriented_faces,is_permitted,newstyles: boolean;
+*)
+
+ var   u,v,du,dv,u1,u2,v1,v2,r1 : real;
+       n1,n2,i,k,ik,iachs,inz,ianf,iplot,nf1,ne0 : integer;
+       ps1,ps2			: vt3d;     {fuer Schnittkanten der Flaechen}
+       inters			: boolean;
+{****************}
+
+{********************************************************************}
+{****  Hilfsprogramme zur Erzeugung von OFF-Dateien fuer GEOMVIEW ***}
+{********************************************************************}
+
+{*****}
+var offdatei : text;
+{*****}
+   
+procedure open_offdatei;
+{ fuer OFF-Datei }
+ var  datei : string;
+ begin
+   writeln('Name der OFF-Datei ? (... .off)');
+   readln(datei);
+   assign(offdatei,datei);
+   rewrite(offdatei);
+ end; { open_offdatei }
+{*************}
+
+procedure close_offdatei;
+ begin
+   Close(offdatei);
+ end;
+{************}
+
+procedure write_nangles_to_offfile;
+var i,k : integer;
+begin
+ open_offdatei;
+ writeln(offdatei, 'OFF');
+ writeln(offdatei,np,' ',nf,' ',0);
+ for i:= 1 to np do with p[i] do
+           writeln(offdatei,'   ',x:3:5,' ',y:3:5,' ',z:3:5);
+ for i:= 1 to nf do
+    with face[i] do
+         begin
+         write(offdatei,npf,'  ');
+         for k:= 1 to npf do write(offdatei,fp[k]-1,' ');
+         writeln(offdatei,' ');
+	 end;   
+ close_offdatei;
+end;  { write_nangles_to_offdatei }
+{*********************************************************************************}
+
+begin {Hauptprogramm}
+   writeln('PLD-Datei ? (ja: 1)');  readln(iplot);
+   graph_on(iplot);
+   repeat
+     writeln('***  Projektion parametrisierter Flaechen         ***');    writeln;
+     writeln('***  hier: "Affensattel" und Zylinder             ***');
+      ne0:= 0;
+{--------------------------------}
+{ Affensattel: Typ ist "quadrangle"}
+     n1:= 30;    n2:= 30;   {Anzahl der Unterteilungen im Parameterbereich};
+     u1:= -0.9;  u2:= 0.9;  {Parametergrenzen}
+     v1:= -0.9;  v2:= 0.9;
+     du:= (u2-u1)/(n1-1);     dv:= (v2-v1)/(n2-1);
+     v:= v1;
+     for k:= 1 to n2 do
+         begin
+         u:= u1;
+         for i:= 1 to n1 do
+             begin
+             ik:=i + (k-1)*n1;
+             put3d(u,v,u*u*u-3*u*v*v, p[ik]);    {Parameterdarstellung}
+             u:= u+du;
+             end;
+         v:= v+dv;
+         end;
+     np:= n1*n2;
+{     aux_quadrangle_triang(n1,n2,true);}
+    aux_quadrangle(n1,n2,0,0,0); 
+     writeln(' np:',np);     writeln(' nf:',nf);     writeln(' ne:',ne);
+     nf1:= nf;
+     for i:= ne0+1 to ne do
+	  with edge[i] do begin color:=cyan; linew:=1; end;
+      ne0:= ne;
+{--------------------------------}      
+{Zylinder 1: x**2 + z**2 -r1**2=0, Typ ist "cylinder"}
+     r1:= 0.6;
+     n1:= 50;    n2:= 10;
+     v1:= -1.5;  v2:= 1;
+     u1:= 0;     u2:= pi2;
+     du:= (u2-u1)/n1;     dv:= (v2-v1)/(n2-1);
+     v:= v1;
+     for k:= 1 to n2 do
+         begin
+         u:= u1;
+         for i:= 1 to n1 do
+             begin
+             ik:=i + (k-1)*n1;
+             put3d(r1*cos(u),r1*sin(u),v, p[np+ik]);
+             u:= u+du;
+             end;
+         v:= v+dv;
+         end;
+     aux_cylinder(n1,n2,np,ne,nf);
+     for i:= ne0+1 to ne do
+	  with edge[i] do begin color:=blue; linew:=1; end;
+      ne0:= ne;
+{---------------------------------------------}
+
+{Durchdringung der beiden Flaechen:}
+    boxes_of_faces;
+    for i:= 1 to nf1 do         { 1.Flaechenschleife }
+        begin
+        for k:= nf1+1 to nf do  { 2.Flaechenschleife }
+            begin
+            is_face_face(i,k, ps1,ps2,inters);
+            if inters then
+               begin
+               p[np+1]:= ps1;  p[np+2]:= ps2;
+               with edge[ne+1] do 
+                    begin ep1:= np+1; ep2:= np+2; vis:= true; end;
+               np:= np+2;  ne:= ne+1;
+               with face[i] do begin fe[nef+1]:= ne; nef:= nef+1;  end;
+               with face[k] do begin fe[nef+1]:= ne; nef:= nef+1;  end;
+               end;  { if }
+            end;  { for k }
+        end;  { for i }
+     for i:= ne0+1 to ne do
+	  with edge[i] do begin color:=red; linew:=1; end;
+      ne0:= ne;
+{---------------------------------------------}
+      
+     repeat
+       writeln;
+       init_centralparallel_projection(2);
+       writeln(' Koordinaten-Achsen ? (Ja = 1)');
+       readln(iachs);
+{ Zeichnen : }
+       draw_area(200,200,100,100,50);
+       if iachs=1 then begin  cp_axes(1.8);  point2d(null2d,0);  end;
+
+       oriented_faces:= false;
+       is_permitted:= true;
+	newstyles:= true;
+       cp_lines_before_convex_faces(oriented_faces,is_permitted,newstyles);
+
+       draw_end ;        writeln ;
+       writeln(' Noch eine Projektion ? ( Ja = 1 )');
+     readln(inz);
+     until inz=0;
+
+{erzeugt off-file fuer GEOMVIEW:}
+   write_nangles_to_offfile;
+
+   writeln('Noch einmal mit ANDEREN  n1, n2 ?  (ja: 1)');
+   readln(ianf);
+   until ianf=0;
+   graph_off;
+ end.
+
+
+
+
+

cdg0gv-2010/beispiele/include

@@ -0,0 +1 @@
+../include

cdg0gv-2010/beispiele/n_eck.p

@@ -0,0 +1,36 @@
+{******************************}
+{***  Regelmaessiges n-Eck  ***}
+{******************************}
+program n_eck;
+uses geograph;
+var    p : vts2d;
+       n,iverb,i,j,inz : integer;
+       r,dw,cdw,sdw : real;
+{*******************}
+ begin {Hauptprogramm}
+   graph_on(0);
+   repeat
+     writeln('***  n-Eck  ***');
+     writeln(' n ?  Radius r des zugehoerigen Kreises ?');
+     readln(n,r);
+     writeln('Jeden Punkt mit jedem Punkt verbinden ? (Ja = 1)');
+     readln(iverb);
+{Berechnung der Eckpunkte:}
+     put2d(r,0, p[0]);
+     dw:= pi2/n;   cdw:= cos(dw);  sdw:= sin(dw);
+     for i := 0 to n-1 do  rotor2d(cdw,sdw, p[i], p[i+1]);
+     draw_area(2*r+20,2*r+20,r+10,r+10,1);
+{Zeichnen:}   new_color(green);
+     if iverb=1 then
+          for i:= 0 to n-1 do
+                for j:= i+1 to n do
+                    line2d(p[i],p[j],0)
+        else
+          curve2d(p,0,n,0);
+     draw_end;
+     writeln('Noch eine Zeichnung? (ja:1, nein:0)');
+     readln(inz);
+   until inz=0;
+   graph_off;
+ end.
+

cdg0gv-2010/beispiele/n_ecki.p

@@ -0,0 +1,41 @@
+{******************************}
+{***  Regelmaessiges n-Eck  ***}
+{******************************}
+program n_eck;
+uses linux;
+const  {$i include/geoconst.pas}
+type   {$i include/geotype.pas}
+var    {$i include/geovar.pas}
+       p : vts2d;
+       n,iverb,i,j,inz : integer;
+       r,dw,cdw,sdw : real;
+     {$i include/geoproc.pas}
+     {$i include/proc_ag.pas}
+{*******************}
+ begin {Hauptprogramm}
+   graph_on(0);
+   repeat
+     writeln('***  n-Eck  ***');
+     writeln(' n ?  Radius r des zugehoerigen Kreises ?');
+     readln(n,r);
+     writeln('Jeden Punkt mit jedem Punkt verbinden ? (Ja = 1)');
+     readln(iverb);
+{Berechnung der Eckpunkte:}
+     put2d(r,0, p[0]);
+     dw:= pi2/n;   cdw:= cos(dw);  sdw:= sin(dw);
+     for i := 0 to n-1 do  rotor2d(cdw,sdw, p[i], p[i+1]);
+     draw_area(2*r+20,2*r+20,r+10,r+10,1);
+{Zeichnen:}   new_color(green);
+     if iverb=1 then
+          for i:= 0 to n-1 do
+                for j:= i+1 to n do
+                    line2d(p[i],p[j],0)
+        else
+          curve2d(p,0,n,0);
+     draw_end;
+    writeln('Noch eine Zeichnung? (ja:1, nein:0)');
+    readln(inz);
+   until inz=0;
+   graph_off;
+ end.
+

cdg0gv-2010/beispiele/pldv

@@ -0,0 +1 @@
+../tools/pldv

cdg0gv-2010/beispiele/povcommandline.txt

@@ -0,0 +1,7 @@
+
+%%%% Type the following command 
+%%%% (file triblend.pp has to be provided by triblend_pov !).
+
+ povray +L/usr/lib/povray31/include +Itriblend.pov +Otriblend.tga +W600 +H400 +FT
+
+%%%% File triblend.tga can be displayed by xv. 

cdg0gv-2010/beispiele/tori.p

@@ -0,0 +1,87 @@
+{***  Projektion zweier Tori mit UP cp_lines_before_convex_faces  ***}
+{********************************************************************}
+program tori_h;
+uses geograph,dos;
+ const {Achtung: es muss  array_size>=nfmax  sein  !!!}
+       nfmax= 20000;  nemax=40000;    nsegmax=10;  npfmax=4;
+
+ 
+ type  vts2d_pol  = array[0..npfmax] of vt2d;
+       vts3d_pol  = array[0..npfmax] of vt3d;
+       r_array_seg  = array[0..nsegmax] of real;
+       i_array_seg  = array[0..nsegmax] of integer;
+       box3d_dat = record
+                   xmin,xmax,ymin,ymax,zmin,zmax : real;
+                   end;
+       face_dat = record
+                  npf,nef : integer;
+                  fp,fe : array[1..npfmax] of integer;
+                  vis : boolean;
+                  box : box3d_dat;
+                  discentre,d : real;
+                  nv : vt3d;
+                  end;
+       edge_dat = record
+                  vis : boolean;
+                  ep1,ep2,color,linew : integer;
+                  end;
+ var   face : array[1..nfmax] of face_dat;
+       edge : array[1..nemax] of edge_dat;
+       p : vts3d;     p0 : vt3d;   pdist : r_array;   error : boolean;
+       n1,n2,np,nf,ne,i,j,k,ik,iachs,inz,ianf,i2tor : integer;
+       r1,r2,xf,yf,x1,y1,x2,y2,dw,cdw,sdw : real;
+     {$i include/proc_zpo.pas}
+
+{****************}
+begin {Hauptprogramm}
+   graph_on(0);
+   repeat
+     writeln(' ***  2 Tori  *** ');    writeln;
+{     writeln(' n2 (>2 Unterteilungen des grossen Kreises) ? ');  readln(n2);
+     writeln(' n1 (>2 Unterteilungen des kleinen Kreises) ? ');  readln(n1);
+     writeln(' 2 Tori ? (ja=1) oder 1 Torus ? ');             readln(i2tor);}
+      n1:= 20; n2:= 20;  i2tor:= 1;
+      r1:= 50;  r2:= 15;
+{ Koordinaten der Punkte 1...n1 (kleiner Kreis): }
+     dw:= pi2/n1;   cdw:= cos(dw);    sdw:= sin(dw);
+     put3d(r1+r2,0,0, p[1]);   put3d(r1,0,0, p0);
+     for i:= 2 to n1 do rotp0y(cdw,-sdw,p0, p[i-1], p[i]);
+{ Koordinaten der restlichen Punkte des 1. Torus: }
+     dw:= pi2/n2;   cdw:= cos(dw);    sdw:= sin(dw);
+     for k:= 2 to n2 do
+         for i:= 1 to n1 do
+             begin
+             ik:= i + (k-1)*n1;
+             rotorz(cdw,sdw,p[ik-n1], p[ik]);
+             end;
+     np:= n1*n2;
+     aux_torus(n1,n2,0,0,0);
+{ 2. Torus:}
+    if i2tor=1 then
+       begin
+       for i:= 1 to np do put3d(p[i].x-r1, -p[i].z, p[i].y,  p[np+i]);
+       aux_torus(n1,n2,np,ne,nf);
+       end;
+     writeln(' np:',np);     writeln(' nf:',nf);     writeln(' ne:',ne);
+
+     repeat
+       init_centralparallel_projection(2);
+      { writeln(' Koordinaten-Achsen ? (Ja = 1)');         readln(iachs);}
+{ Zeichnen : }
+	draw_area(200,200,80,90,1);
+    {   if iachs=1 then begin  cp_axes(20);  point2d(null2d,0);  end; }
+
+       cp_lines_before_convex_faces(true,false,false);    {Hiddenline-Algorithmus}
+
+       draw_end ;        writeln ;
+	writeln(' Noch eine Projektion ? ( Ja = 1 )');     readln(inz);
+     until inz=0;
+   writeln('Noch einmal mit ANDEREN  n1, n2 ?  (ja: 1)'); readln(ianf);
+   until ianf=0;
+   graph_off;
+ end.
+
+
+
+
+

cdg0gv-2010/beispiele/triang_proc.p

@@ -0,0 +1,718 @@
+
+(*
+
+{***************************************************}
+{*********  P R O C E D U R E S  ON VECTORS  *******}
+{***************************************************}
+
+const array_size = 20000; 
+      pi= 3.14159265358;       pi2= 6.2831853;       pih= 1.5707963;
+      eps8=0.00000001;
+
+type  vt3d    = record  x: real;  y: real;  z: real; end;
+      vts3d   = array[0..array_size] of vt3d;
+{***************************************************}
+
+function max(a,b : real) : real;
+ begin   if a>=b then  max:= a   else   max:= b;   end; {max}
+{*************}
+function min(a,b : real) : real;
+ begin   if a<=b then  min:= a   else   min:= b;   end; {min}
+{*************}
+
+procedure put3d(x,y,z : real; var v: vt3d);
+ begin   v.x:= x;  v.y:= y;  v.z:= z;   end;
+{*************}
+
+procedure get3d(v : vt3d; var x,y,z: real);
+ begin  x:= v.x;  y:= v.y;  z:= v.z;  end;
+{*************}
+
+procedure scale3d(r : real; v: vt3d; var vs: vt3d);
+ begin   vs.x:= r*v.x;  vs.y:= r*v.y;  vs.z:= r*v.z;  end;
+{*************}
+
+procedure sum3d(v1,v2 : vt3d; var vs : vt3d);
+ begin vs.x:= v1.x + v2.x;  vs.y:= v1.y + v2.y;  vs.z:= v1.z + v2.z; end;
+{*************}
+
+procedure diff3d(v1,v2 : vt3d; var vs : vt3d);
+ begin vs.x:= v1.x - v2.x;  vs.y:= v1.y - v2.y;  vs.z:= v1.z - v2.z; end;
+{*************}
+
+procedure lcomb2vt3d(r1: real; v1: vt3d; r2: real; v2: vt3d; var vlc : vt3d);
+ begin
+ vlc.x:= r1*v1.x + r2*v2.x; vlc.y:= r1*v1.y + r2*v2.y; vlc.z:= r1*v1.z + r2*v2.z;
+ end;
+{*************}
+
+procedure lcomb3vt3d(r1: real; v1: vt3d; r2: real; v2: vt3d;
+                                          r3: real; v3: vt3d; var vlc : vt3d);
+ begin
+   vlc.x:= r1*v1.x + r2*v2.x + r3*v3.x;
+   vlc.y:= r1*v1.y + r2*v2.y + r3*v3.y;
+   vlc.z:= r1*v1.z + r2*v2.z + r3*v3.z;
+ end;
+{*************}
+
+function abs3d(v : vt3d) : real;
+ begin  abs3d:= abs(v.x) + abs(v.y) + abs(v.z);  end;
+{*************}
+function length3d(v : vt3d) : real;
+ begin  length3d:= sqrt( sqr(v.x) + sqr(v.y) + sqr(v.z));  end;
+{*************}
+procedure normalize3d(var p: vt3d);
+var c : real;
+begin c:= 1/length3d(p); p.x:= c*p.x; p.y:= c*p.y; p.z:= c*p.z  end;
+{************}
+function scalarp3d(p1,p2 : vt3d) : real;
+ begin  scalarp3d:= p1.x*p2.x + p1.y*p2.y + p1.z*p2.z;  end;
+{*************}
+function distance3d(p,q : vt3d): real;
+begin  distance3d:= sqrt( sqr(p.x-q.x) + sqr(p.y-q.y) + sqr(p.z-q.z) );  end;
+{*************}
+function distance3d_square(p,q : vt3d) : real;
+begin distance3d_square:= sqr(p.x-q.x) + sqr(p.y-q.y) + sqr(p.z-q.z);  end;
+{**************}
+
+procedure vectorp(v1,v2 : vt3d; var vp : vt3d);
+{vectorproduct of v1,v2}
+ begin
+   vp.x:=  v1.y*v2.z - v1.z*v2.y ;
+   vp.y:= -v1.x*v2.z + v2.x*v1.z ;
+   vp.z:=  v1.x*v2.y - v2.x*v1.y ;
+ end; {vectorp}
+{*************}
+
+function determ3d(v1,v2,v3: vt3d) : real;
+{determinant of 3 vectors.}
+ begin
+   determ3d:= v1.x*v2.y*v3.z + v1.y*v2.z*v3.x + v1.z*v2.x*v3.y
+            - v1.z*v2.y*v3.x - v1.x*v2.z*v3.y - v1.y*v2.x*v3.z;
+ end; {determ3d}
+{*************}
+
+procedure rotorz(cos_rota,sin_rota : real; p : vt3d; var pr: vt3d);
+ begin
+   pr.x:= p.x*cos_rota - p.y*sin_rota;
+   pr.y:= p.x*sin_rota + p.y*cos_rota;    pr.z:= p.z;
+ end; {rotorz}
+{*************}
+
+function polar_angle(x,y : real) : real;
+ {determines the polar angle of point (x,y)}
+ var w : real;
+ begin
+   if (x=0) and (y=0) then   w:= 0
+     else
+      begin
+        if abs(y)<=abs(x) then
+           begin
+             w:= arctan(y/x);
+             if x<0 then  w:=pi+w
+                else  if (y<0) and( w<>0) then  w:= pi2+w;
+           end
+          else
+           begin
+             w:= pih-arctan(x/y);
+             if y<0 then w:= pi+w;
+           end; {if}
+      end;  {if}
+    polar_angle:= w;
+ end;  { polar_angle }
+{******************}
+
+procedure newcoordinates3d(p,b0,b1,b2,b3: vt3d; var pnew: vt3d);
+{determines the coordinates of p for the basis b1,b2,b3 with origin b0.}
+var det : real;   p0: vt3d;
+begin
+  diff3d(p,b0, p0);              det:= determ3d(b1,b2,b3);
+  pnew.x:= determ3d(p0,b2,b3)/det;
+  pnew.y:= determ3d(b1,p0,b3)/det;
+  pnew.z:= determ3d(b1,b2,p0)/det;
+end;  { newcoordinates3d }  
+{************}
+
+*)
+
+{*******************************************************************************}
+
+{**********************************************}
+{********** T R I A N G U L A T I O N : *******}     
+{**********************************************}
+
+(*
+{for triangulation:}
+const tnfmax=10000; tnemax=20000; tnpmax=10000; 
+      tfrnpmax= 500;  {max. number of points in a front}  
+      tnfrmax=20;     {max. number of further fronts.}*)
+
+type
+(*      implicit3d = procedure(p: vt3d; var fvalue: real; var gradf: vt3d);*)
+
+      tface_dat = record
+                  p1,p2,p3: integer;        {points of a triangle}
+                  end;
+      tpoint_dat= record
+                  p,nv,tv1,tv2 : vt3d;      {coordinates, normal,tangentvectors}
+		     full,achange: boolean; {full=true: point surrounded by triang.} 
+		     angle: real;           {achange=true: angle was changed}
+                  end;
+var
+      tface : array[1..tnfmax] of tface_dat;
+      tpoint: array[1..tnpmax] of tpoint_dat;           
+      tfrontpt: array[1..tfrnpmax] of integer; {points in actual front polygon}  
+      tfrontnp: integer;     {number of front points} 
+      tfr : array[0..tnfrmax,0..tfrnpmax] of integer;   {further fronts} 
+      tfrbox : array[0..tnfrmax,1..6] of real; 
+      tnp,tne,tnf,tnfr : integer;  {number of points, edges, triangles, fronts}
+      tstepl,tstepl_square,xmin,xmax,ymin,ymax,zmin,zmax,minangle : real;
+      n_triang,fullcount,frontnumber: integer;
+      f_gradf : implicit3d;
+      actual_ip : integer;
+      nearpointtest : boolean;
+      cuttype : integer;
+      xcut,ycut,zcut,rcut,rcut_square: real;
+{*****************************************}
+procedure quader(xmin,xmax,ymin,ymax,zmin,zmax: real);
+var p1,p2,p3,p4,p5,p6,p7,p8: vt3d;
+begin
+  put3d(xmin,ymin,zmin, p1);
+  put3d(xmax,ymin,zmin, p2);
+  put3d(xmax,ymax,zmin, p3);
+  put3d(xmin,ymax,zmin, p4);
+  put3d(xmin,ymin,zmax, p5);
+  put3d(xmax,ymin,zmax, p6);
+  put3d(xmax,ymax,zmax, p7);
+  put3d(xmin,ymax,zmax, p8);
+  cp_line(p1,p2,0);  cp_line(p2,p3,0);  cp_line(p3,p4,0);  cp_line(p4,p1,0);
+  cp_line(p5,p6,0);  cp_line(p6,p7,0);  cp_line(p7,p8,0);  cp_line(p8,p5,0);
+  cp_line(p1,p5,0);  cp_line(p2,p6,0);  cp_line(p3,p7,0);  cp_line(p4,p8,0);
+end;  
+{***************************************************************************}
+
+
+procedure insert_point(ip_new,i_new: integer);
+{inserts point ip_new as i_new-th front point}
+var i : integer;
+begin
+  for i:= tfrontnp+1 downto i_new do tfrontpt[i+1]:=  tfrontpt[i];
+  tfrontpt[i_new]:= ip_new; 
+  tfrontnp:= tfrontnp+1;
+end;
+{******} 
+procedure delete_point(i_del: integer);
+{deletes i_del-th front point}
+var i : integer;
+begin
+  tfrontnp:= tfrontnp-1;
+  for i:= i_del to tfrontnp+1 do tfrontpt[i]:=  tfrontpt[i+1];
+end;
+{******} 
+
+function  point_ok_box(p: vt3d): boolean;
+{point p in bounding box ?}
+begin
+   with p do
+      if ((x<=xmax) and (x>=xmin) and 
+	  (y<=ymax) and (y>=ymin) and 
+	  (z<=zmax) and (z>=zmin)) then point_ok_box:= true
+	                           else point_ok_box:= false;
+end;
+{*****}
+procedure cut_seg_box(p0,p: vt3d; var pc: vt3d);
+{cuts new edge p0-p at bounding box -> p0-pc}
+var dv: vt3d; x0,y0,z0,t: real;
+begin
+   diff3d(p,p0, dv);  t:= 1;  get3d(p0, x0,y0,z0);
+   with dv do
+          begin 
+	     if x>0 then t:= min(t,(xmax-x0)/x) else t:= min(t,(xmin-x0)/x);
+	     if y>0 then t:= min(t,(ymax-y0)/y) else t:= min(t,(ymin-y0)/y); 
+	     if z>0 then t:= min(t,(zmax-z0)/z) else t:= min(t,(zmin-z0)/z);
+	  end; 
+ lcomb2vt3d(1,p0,t,dv, pc); 
+end;
+{******}     
+
+function point_ok_cyl(p: vt3d): boolean;
+{point p in bounding box ?}
+begin
+   with p do
+	if (z<=zmax) and (z>=zmin) and (sqr(x-xcut)+sqr(y-ycut)< rcut_square)
+	   then point_ok_cyl:= true
+	   else point_ok_cyl:= false;
+end;
+{*****}
+procedure cut_seg_cyl(p0,p: vt3d; var pc: vt3d);
+{cuts new edge p0-p at bounding cylinder -> p0-pc}
+var dv: vt3d; t,x0,y0,z0,a,b,d2: real;
+begin
+   diff3d(p,p0, dv);  t:= 1;  get3d(p0, x0,y0,z0);
+   with dv do
+          begin 
+	  if z>0 then t:= min(t,(zmax-z0)/z) else t:= min(t,(zmin-z0)/z);
+          d2:=x*x+y*y;
+          if d2>0 then
+             begin
+             x0:= p0.x-xcut;  y0:=p0.y-ycut;   
+             a:= (x0*x+y0*y)/d2;  b:= (x0*x0+y0*y0-rcut_square)/d2;
+             t:= min(t,-a+sqrt(a*a-b)); 
+             end;
+	  end; 
+   lcomb2vt3d(1,p0,t,dv, pc); 
+end;
+{******}  
+function point_ok_sph(p: vt3d): boolean;
+{point p in bounding sphere ?}
+begin
+   with p do
+	if (sqr(x-xcut)+sqr(y-ycut)+sqr(z-zcut)< rcut_square)
+	   then point_ok_sph:= true
+	   else point_ok_sph:= false;
+end;
+{*****}
+procedure cut_seg_sph(p0,p: vt3d; var pc: vt3d);
+{cuts new edge p0-p at bounding sphere -> p0-pc}
+var dv: vt3d; t,x0,y0,z0,a,b,d2: real;
+begin
+   diff3d(p,p0, dv);  t:= 1;  get3d(p0, x0,y0,z0);
+   with dv do
+          begin 
+          d2:=x*x+y*y+z*z;
+          if d2>0 then
+             begin
+	     x0:= p0.x-xcut;  y0:=p0.y-ycut; z0:= p0.z-zcut;
+             a:= (x0*x+y0*y+z0*z)/d2;
+	     b:= (x0*x0+y0*y0+z0*z0-rcut_square)/d2;
+             t:= min(t,-a+sqrt(a*a-b)); 
+             end;
+	  end; 
+   lcomb2vt3d(1,p0,t,dv, pc); 
+end;
+{******}  
+
+function point_ok(p: vt3d): boolean;
+{point p in bounding box or cylinder}
+begin
+  if cuttype=1 then
+     point_ok:= point_ok_cyl(p)
+  else if cuttype=2 then point_ok:= point_ok_sph(p)
+                    else point_ok:= point_ok_box(p);
+end;
+{*****}
+procedure cut_seg(p0,p: vt3d; var pc: vt3d);
+{cuts new edge p0-p at bounding box or cylinder -> p0-pc}
+begin
+  if cuttype=1 then cut_seg_cyl(p0,p, pc)
+  else if cuttype=2 then cut_seg_sph(p0,p, pc)
+                    else cut_seg_box(p0,p, pc);
+end;
+{******}  
+
+procedure surface_point_normal_tangentvts(p_start : vt3d;
+                                          f_gradf : implicit3d;  
+                                  var p,nv,tv1,tv2: vt3d);
+{***}
+procedure surface_point(p0 : vt3d; var p_surface,nv : vt3d);
+{seeks surface point along the steepest way}
+ var delta,fi : real;   p_i,p_i1,dv,gradfi : vt3d;  
+{**}
+procedure newton_step(p_i : vt3d; var p_i1 : vt3d);
+ var  t,cc : real;  
+ begin
+   f_gradf(p_i, fi,gradfi);
+   cc:= scalarp3d(gradfi,gradfi);
+    if cc>1E-15{eps8} then t:= -fi/cc
+    else
+    begin t:= 0; writeln(cc,' WARNING tri (surface_point...): newton');
+       new_color(lightred);  cp_point(p_i, 1); new_color(default); end;
+   lcomb2vt3d(1,p_i, t,gradfi, p_i1);
+ end;
+{**}
+ begin {surface_point}
+   p_i:= p0;      
+   repeat
+     newton_step(p_i, p_i1);     diff3d(p_i1,p_i, dv);
+     delta:= abs3d(dv);
+     p_i:= p_i1;
+   until delta < eps8;
+   p_surface:= p_i1; nv:= gradfi;  
+ end;   { surface_point }
+{***}
+begin
+  surface_point(p_start, p,nv);
+  normalize3d(nv);
+  with nv do if (abs(x)>0.5) or (abs(y)>0.5) then put3d(y,-x,0, tv1)
+                                             else put3d(-z,0,x, tv1);
+  normalize3d(tv1);  vectorp(nv,tv1, tv2);  {nv,tv1,tv2: ON-base} 
+end;  {surface_point_normal_tangentvts}
+{****}
+
+procedure start_triangulation(p_start: vt3d; f_gradf: implicit3d);
+{calculates from p_start the first surface point p[1] and the first six triangles}
+var p0,p1,nv0,tv10,tv20 : vt3d;    i,tnp0 : integer;   dw,cw,sw : real;
+begin
+  tnp:= tnp+1; 
+  surface_point_normal_tangentvts(p_start,f_gradf,p0,nv0,tv10,tv20);
+  if not point_ok(p0) then
+  begin writeln('!!! First point not in bounding box or cylinder !!!');
+	n_triang:= 0;  exit;
+     end; 
+  with tpoint[tnp] do 
+  begin p:= p0;  nv:= nv0;  tv1:= tv10;  tv2:= tv20;  end;
+  tnp0:= tnp;
+  for i:= 0 to 5 do
+      begin
+      dw:= pi/3; cw:= cos(i*dw);  sw:= sin(i*dw);
+      lcomb3vt3d(1,p0,tstepl*cw,tv10,tstepl*sw,tv20, p1);
+      tnp:= tnp +1;  tfrontpt[i+1]:= tnp;  
+      with tpoint[tnp] do
+           begin 
+           surface_point_normal_tangentvts(p1,f_gradf, p,nv,tv1,tv2);
+           achange:= true;
+           end;
+      end;   { for }
+  tfrontnp:= 6;   tpoint[1].full:= true;
+  for i:= 1 to 6 do   {triangles}
+      begin
+      tnf:= tnf+1;                 
+      with tface[tnf] do begin p1:= tnp0; p2:= tnp0+i; p3:= tnp0+i+1; end;
+      end; 
+  tface[tnf].p3:= tnp0+1;       
+end;  
+{********}  
+procedure new_triangle(q1,q2,q3: integer);
+begin
+  tnf:= tnf+1;
+  with tface[tnf] do begin p1:= q1; p2:= q2; p3:= q3; end;
+end;   
+{****}
+function reduce(n: integer): integer;
+begin
+  reduce:= n;
+  if n<1        then reduce:= n+tfrontnp;
+  if n>tfrontnp then reduce:= n-tfrontnp;
+end;
+{****}
+
+procedure make_angle(ipf: integer);
+{calculates the front angle at ipf-th front point}
+var pn1,pn2,pn11,pn22 : vt3d;    ip,ip1,ip2 : integer;    w1,w2: real; 
+begin
+  ip:=  tfrontpt[ipf];    
+  with tpoint[ip] do
+       if (full or (not achange) ) then exit
+          else
+           begin
+	   ip1:= tfrontpt[reduce(ipf-1)];  
+           ip2:= tfrontpt[reduce(ipf+1)];  
+           pn1:= tpoint[ip1].p;  newcoordinates3d(pn1,p,tv1,tv2,nv, pn11);  
+           pn2:= tpoint[ip2].p;  newcoordinates3d(pn2,p,tv1,tv2,nv, pn22); 
+           w1:= polar_angle(pn11.x,pn11.y);    
+           w2:= polar_angle(pn22.x,pn22.y);     if w2<w1 then w2:= w2+pi2;
+           angle:= w2-w1;   achange:= false;
+           end;
+end;  
+{****}
+
+function fangle(ipf: integer): real;
+{calculates the front angle at ipf-th front point}
+var pn1,pn2,pn11,pn22 : vt3d;    ip,ip1,ip2 : integer;    w1,w2: real; 
+begin
+  ip:=  tfrontpt[ipf];    
+  with tpoint[ip] do
+       begin
+       ip1:= tfrontpt[reduce(ipf-1)]; 
+       ip2:= tfrontpt[reduce(ipf+1)];  
+       pn1:= tpoint[ip1].p;  newcoordinates3d(pn1,p,tv1,tv2,nv, pn11);  
+       pn2:= tpoint[ip2].p;  newcoordinates3d(pn2,p,tv1,tv2,nv, pn22); 
+       w1:= polar_angle(pn11.x,pn11.y);    
+       w2:= polar_angle(pn22.x,pn22.y);     if w2<w1 then w2:= w2+pi2;
+       fangle:= w2-w1; 
+       end;
+end;  
+{****}
+
+function outside(ip0,ip1,iptest: integer): boolean;
+{is point p[iptest] in not triangulated region at point p[ip0] ?}
+var pn1,pn11,ptest,ptestt: vt3d;  w1,wtest: real; 
+begin
+  with tpoint[ip0] do
+       begin
+       pn1:= tpoint[ip1].p;       newcoordinates3d(pn1,p,tv1,tv2,nv, pn11);  
+       ptest:= tpoint[iptest].p;  newcoordinates3d(ptest,p,tv1,tv2,nv, ptestt); 
+       w1:= polar_angle(pn11.x,pn11.y);    
+       wtest:= polar_angle(ptestt.x,ptestt.y);     
+       if wtest<w1 then wtest:= wtest+pi2;
+       if wtest<w1+angle then outside:= true else outside:= false; 
+       end;
+end;   
+{********}  
+
+procedure test_criticalpt(ipf: integer; var ipf_nearp: integer);
+{seeks another front point near to ipf-th front point}
+label 10;
+var i,k,ip,ip1,ip2,ip11,ip22,ip111,ip222,ipi,nn: integer;    p0: vt3d;
+begin
+  ip:= tfrontpt[ipf];                ipf_nearp:= -1;    
+  with tpoint[ip] do if full then exit else p0:= p; 
+  if ipf<3 then nn:= tfrontnp-3+ipf else nn:= tfrontnp;
+  {tests only front points i>ipf+3 ! important for divide_front}     
+  for i:= ipf+3 to nn do
+      begin               {ipi not neighbor or neighbor of neighbors} 
+      ipi:= tfrontpt[i];     
+      with tpoint[ipi] do 
+           if not full then 
+              if abs(p0.x-p.x)<tstepl then
+                 if abs(p0.y-p.y)<tstepl then
+                    if abs(p0.z-p.z)<tstepl then    
+                       if distance3d_square(p0,p) < tstepl_square then 
+                          begin 
+                          if not outside(ip,tfrontpt[reduce(ipf-1)],ipi) 
+                                                               then goto 10;
+        		  if scalarp3d(tpoint[ip].nv,tpoint[ipi].nv)<0
+			                                       then goto 10;
+                          ipf_nearp:= i;   exit; 
+                          end; { if }
+10:   end;  { for }          
+end;  { test_criticalpt }
+{*********}
+
+procedure find_pair_of_nearpts(var ipf,ipf_nearp,frontnumber: integer);
+{seeks a nearpoint to a front point}
+label 10;
+var i1,i2,i,k,j,l,ip,ipj,inp,tfrontnpi: integer;   p0: vt3d; 
+begin
+  ipf:= 0;  ipf_nearp:= -1;  frontnumber:= 0;
+  for i:= 1 to tfrontnp  do 
+        begin 
+        test_criticalpt(i,inp);  
+        if inp>0 then 
+           begin ipf:= i; ipf_nearp:= inp; frontnumber:= 0; exit; end;  
+        end;
+  for i:= 1 to tnfr do
+      begin
+      tfrontnpi:= tfr[i,0];       
+      for k:= 1 to tfrontnp do
+          begin
+          ip:=  tfrontpt[k];  
+          with tpoint[ip] do if full then goto 10 else p0:= p;   
+          for j:= 1 to tfrontnpi do
+              begin 
+              ipj:= tfr[i,j];     
+              with tpoint[ipj] do 
+                   if abs(p0.x-p.x)<tstepl then
+                      if abs(p0.y-p.y)<tstepl then
+                         if abs(p0.z-p.z)<tstepl then 
+                            if distance3d_square(p0,p) < tstepl_square then 
+                               if outside(ip,tfrontpt[reduce(k-1)],ipj) then 
+                                   begin
+                                   frontnumber:= i;  ipf:= k;  ipf_nearp:= j;       
+                                   exit; 
+                                   end;
+              end;  { for j }                 
+10:	  end;  { for k }
+      end;   { for i }                
+end;                    
+{**************}
+
+function minanglept: integer;
+{determines the front point with minimal front angle.}
+var  i,imin: integer; min: real;
+begin
+  min:= 10; 
+  for i:= 1 to tfrontnp do with tpoint[tfrontpt[i]] do
+      if ((not full) and (angle<min)) 
+         then begin min:= angle; imin:= i; end;
+  if min<10 then begin minanglept:= imin; minangle:= min; end 
+            else minanglept:= -1;      
+end;
+{****}
+
+procedure complete_point(ipf: integer; f_gradf: implicit3d);
+{determines for ipf-th front point the missing triangles}
+var p0,nv0,tv10,tv20,pn1,pn2,pn0,p_start,pc_start,pn11,pn22 : vt3d;    
+    i,ip,ip1,ip2,ipf1,ipf2,ne_rest : integer;      dw,cdw,sdw: real; 
+begin
+  if tfrontnp<=3 then exit;
+  ip:= tfrontpt[ipf];
+  actual_ip:= ip;   {global var. for start parameters of footpoint algorithms}
+  with tpoint[ip] do
+       begin 
+       if full then exit;
+       ipf1:= reduce(ipf-1);  ipf2:= reduce(ipf+1);
+       ip1:= tfrontpt[ipf1];  ip2:= tfrontpt[ipf2];  
+       pn1:= tpoint[ip1].p;  newcoordinates3d(pn1,p,tv1,tv2,nv, pn11);  
+       pn11.z:=0; normalize3d(pn11);  scale3d(tstepl,pn11, pn11);
+       pn2:= tpoint[ip2].p;  newcoordinates3d(pn2,p,tv1,tv2,nv, pn22); 
+       if achange then angle:=fangle(ipf);
+       ne_rest:= trunc(angle*3/pi);
+       dw:= angle/(ne_rest+1);   
+       if (dw<0.8) and (ne_rest>0) 
+          then begin ne_rest:= ne_rest-1;  dw:= angle/(ne_rest+1); end;
+       if (ne_rest=0) and (dw>0.8) and (distance3d(pn1,pn2)>1.25*tstepl) then 
+           begin ne_rest:= 1;  dw:= dw/2; end;       
+       p0:= p; nv0:= nv; tv10:= tv1; tv20:= tv2; 
+       end;  
+  if ((distance3d_square(p0,pn1)<0.2*tstepl_square) or    
+      (distance3d_square(p0,pn2)<0.2*tstepl_square))  then ne_rest:=0;  
+  if ne_rest=0 
+     then new_triangle(ip1,ip2,ip)
+     else   
+     for i:= 1 to ne_rest do
+         begin
+         cdw:= cos(dw);  sdw:= sin(dw);
+         rotorz(cdw,sdw,pn11, pn11);
+         lcomb3vt3d(1,p0,pn11.x,tv10,pn11.y,tv20, p_start);
+         tnp:= tnp +1;  
+	 if point_ok(p_start)
+	    then begin tpoint[tnp].full:= false; pc_start:=p_start; end
+            else begin tpoint[tnp].full:= true; cut_seg(p0,p_start, pc_start); end;
+         with tpoint[tnp] do
+              begin                           
+              surface_point_normal_tangentvts(pc_start,f_gradf, p,nv,tv1,tv2);
+	      if scalarp3d(nv,tpoint[ip].nv)<0 then  {change if necessary (SING.)}
+                                  begin scale3d(-1,nv,nv); change3d(tv1,tv2); end;
+	      achange:= true;
+              end;
+         if i=1 then new_triangle(ip1,tnp,ip); 
+         if i=ne_rest then new_triangle(tnp,ip2,ip) 
+                      else new_triangle(tnp,tnp+1,ip);
+         end;   { for }
+  delete_point(ipf);      
+  for i:= 0 to ne_rest-1 do insert_point(tnp-i,ipf);
+  tpoint[ip1].achange:= true;  tpoint[ip2].achange:= true;
+end;  
+{********}  
+
+procedure divide_front(ipf1,ipf2: integer);
+var i,nn: integer;   fa1,fa2: real;  
+begin 
+  tnfr:= tnfr+1;
+  with tpoint[tfrontpt[ipf1]] do achange:= true; 
+  with tpoint[tfrontpt[ipf2]] do achange:= true; 
+  for i:= 0 to ipf2-ipf1 do tfr[tnfr,i+1]:= tfrontpt[ipf1+i];    
+  for i:=1 to tfrontnp-ipf2+1 do tfrontpt[ipf1+i]:= tfrontpt[ipf2+i-1];
+  tfrontnp:= tfrontnp-(ipf2-ipf1-1);
+  tfr[tnfr,0]:=  ipf2-ipf1+1;
+  fa1:= fangle(ipf1); fa2:= fangle(ipf1+1);  
+  if fa1<fa2 then begin nn:= tfrontnp; 
+                        complete_point(ipf1  ,f_gradf);  
+                        complete_point(ipf1+tfrontnp-nn+1  ,f_gradf); end
+             else begin complete_point(ipf1+1,f_gradf); 
+                        complete_point(ipf1  ,f_gradf); end;
+end; 
+{***********}
+
+procedure unite_front(ipf1,ipf2,frontnumber: integer);
+{unites the actual front polygon with a further front (frontnumber)}
+var i,ipf,tfrontnpi,nn: integer;   fa1,fa2: real;
+begin
+  tfrontnpi:=  tfr[frontnumber,0];
+  for i:= 1 to tfrontnpi do 
+      begin
+      ipf:= ipf2+i-1;
+      if ipf>tfrontnpi then ipf:= ipf-tfrontnpi; 
+      insert_point(tfr[frontnumber,ipf],ipf1+i);
+      end;
+  insert_point(tfrontpt[ipf1+1],ipf1+tfrontnpi+1);
+  insert_point(tfrontpt[ipf1],ipf1+tfrontnpi+2);
+  for i:= 0 to tfr[tnfr,0] do tfr[frontnumber,i]:= tfr[tnfr,i];
+  tnfr:= tnfr-1;
+  with tpoint[tfrontpt[ipf1]] do 
+       begin achange:= true; full:= false; end;
+  with tpoint[tfrontpt[ipf1+1]] do 
+       begin achange:= true; full:= false; end;
+  fa1:= fangle(ipf1); fa2:= fangle(ipf1+1);  
+  if fa1<fa2 then begin nn:= tfrontnp; 
+                        complete_point(ipf1  ,f_gradf);  
+                        complete_point(ipf1+tfrontnp-nn+1  ,f_gradf); end
+             else begin complete_point(ipf1+1,f_gradf); 
+                        complete_point(ipf1  ,f_gradf); end;
+end;  { unite_front }    
+{*******}
+
+procedure triangulation(nearpointtest : boolean);
+var i,k,ipi1,ipi2,ipf,ipf1,ipf2,delay: integer;
+begin       
+{Box-Tests:}
+  for i:= 1 to tfrontnp do
+      if not point_ok(tpoint[tfrontpt[i]].p) then
+          begin writeln('!!! STARTING front points NOT in bounding BOX !!!');
+          quader(xmin,xmax,ymin,ymax,zmin,zmax);
+	     new_linewidth(10); cp_point(tpoint[tfrontpt[i]].p,0); new_linewidth(1); 
+          n_triang:= 0;  exit;
+          end;
+  for k:= 1 to tnfr do 
+  for i:= 1 to tfr[k,0] do 
+      if not point_ok( tpoint[tfr[k,i]].p) then
+          begin writeln('!!! BOUNDARY points NOT in bounding BOX !!!');
+          quader(xmin,xmax,ymin,ymax,zmin,zmax);
+	  new_linewidth(10); cp_point(tpoint[tfr[k,i]].p,0); new_linewidth(1);
+          n_triang:= 0;  exit;
+          end; 
+  writeln('**********************************'); 
+   if nearpointtest then writeln('triang.: nearpoint test is ON !');
+   if cuttype=1 then writeln('triang.: bounded by CYLINDER ')
+       else if cuttype=2 then writeln('triang.: bounded by SPHERE ')
+                         else writeln('triang.: bounded by a BOX');
+  writeln('**********************************'); 
+  tstepl_square:= tstepl*tstepl; fullcount:= 0;
+  while (tnfr>=0) and (tnf<n_triang) and (fullcount<tfrontnp) do
+        begin 
+	   fullcount:= 0;  delay:= 0;   minangle:= 0;
+        while (tnf<n_triang) and (fullcount<tfrontnp) and (tfrontnp>3)  do   
+	      begin             
+	      if ((not nearpointtest) or (tfrontnp<10) or (minangle<1.5))
+                   then delay:= 1;
+              if delay=0 then find_pair_of_nearpts(ipf1,ipf2,frontnumber)
+		         else begin delay:= delay-1; ipf2:= -1; end;		   
+              if ipf2>0 then
+                   begin
+                   ipi1:= tfrontpt[ipf1];
+                   if frontnumber=0 then  ipi2:= tfrontpt[ipf2]
+                                    else  ipi2:= tfr[frontnumber,ipf2];
+                   delay:= 0;
+		   writeln('triangles: ',tnf,'  NEARPOINTS detected: ',ipi1,' +++ ',ipi2);
+                   if frontnumber=0 then divide_front(ipf1,ipf2); 
+                   if frontnumber>0 then unite_front(ipf1,ipf2,frontnumber);
+                   end;
+             for i:= 1 to tfrontnp do make_angle(i); 
+             ipf:= minanglept;   
+             if ipf>0 then complete_point(ipf,f_gradf);
+             for i:= 1 to tfrontnp do make_angle(i); 
+	     ipf:= minanglept; {renews minangle for "if ....(minangle<..)" above !!}
+             fullcount:= 0;      
+             for i:= 1 to tfrontnp do  
+                 with tpoint[tfrontpt[i]] do if full then fullcount:= fullcount+1;
+             end;  { while }
+             
+        if tfrontnp=3 then 
+           begin
+           new_triangle(tfrontpt[1],tfrontpt[2],tfrontpt[3]);
+           tfrontnp:= 0;
+           end;
+        if ((tfrontnp=0) or (fullcount=tfrontnp)) and (tnfr>0) then 
+           begin
+           tfrontnp:= tfr[tnfr,0];  
+           for i:=1 to tfrontnp do 
+               begin
+               tfrontpt[i]:= tfr[tnfr,i];
+               with tpoint[tfrontpt[i]] do achange:= true;
+               end;
+           tnfr:= tnfr-1;   fullcount:=0;   writeln('tnfr: ',tnfr);   
+           end;
+	end;  { while }
+  writeln('*****************************************************'); 
+  writeln('triang.: total number of triangels: ', tnf);  
+  writeln('triang.: remaining front points: ', tfrontnp);  
+  writeln('triang.: remaining fronts: ', tnfr);  
+  writeln('*****************************************************'); 
+   if tnf>tnfmax then writeln('triang. warning: tnfmax to small !!!');
+   if tnp>tnpmax then writeln('triang. warning: tnpmax to small !!!');
+   if tfrontnp>tfrnpmax  then writeln('triang. warning: tfrnpmax to small !!!');
+   if tnfr>tnfrmax  then writeln('triang. warning: tnfrmax to small !!!');
+end;   { triangulation }
+{*************************************************************}
+

cdg0gv-2010/beispiele/triblend.pov

@@ -0,0 +1,93 @@
+// Persistence Of Vision raytracer version 3.0 sample file.
+// PolyWood.pov - Wooden polyhedron hollowed by a sphere
+// on a grassy hilly lawn.  Shows how easy it is to
+// create interesting shapes with CSG operations
+// and simple primitive shapes.
+// File by Eduard [esp] Schwan
+
+#version 3.0
+global_settings { assumed_gamma 1.0 }
+
+#include "shapes.inc"
+#include "colors.inc"
+#include "textures.inc"
+
+// Moi
+camera {
+   location   <7,  2.5, -8.0>
+   direction  <0.0,   0.0,   1.0>
+   up         <0.0,   1.0,   0.0>
+   right      <4/3,   0.0,   0.0>
+   look_at    <0,     0,     0>
+}
+
+// Some Light just above the horizon for a long shadow
+light_source
+{
+  <20, 5, -60>
+  color White
+}
+light_source{<60, 5, -20> color White}
+
+
+// The Cloudy Blue Sky
+sphere
+{
+  <0, 0, 0>, 10000
+  pigment
+  {
+     Bright_Blue_Sky
+     scale <4000, 600, 1000>
+  }
+}
+
+
+// The Hilly Grassy Land
+plane
+{
+  y, -3
+  pigment { color red 0.2 green 3.0 blue 0.4 }
+  finish
+  {
+      crand 0.025 // a little randomness to hide the rather severe color banding
+      ambient 0.1
+      diffuse 0.7
+      roughness 1
+   }
+   normal { bumps 0.5  scale 10 }
+}
+
+union{
+       #include "triblend.pp"
+       rotate <-90, 15, 0>
+       translate <1, 0, 0>
+       pigment{color Green}
+   pigment {
+      wood
+      turbulence 0.04
+      colour_map {
+        [0.0 0.4  color red 0.8 green 0.4 blue 0.2
+                  color red 0.8 green 0.4 blue 0.1]
+        [0.4 0.5  color red 0.1 green 0.3 blue 0.1
+                  color red 0.1 green 0.3 blue 0.2]
+        [0.5 0.8  color red 0.1 green 0.3 blue 0.2
+                  color red 0.8 green 0.4 blue 0.1]
+        [0.8 1.0  color red 0.8 green 0.4 blue 0.1
+                  color red 0.8 green 0.4 blue 0.2]
+                  }
+      scale <0.2, 0.2, 1>
+      rotate <45, 0, 5>
+      translate <2, 2, -4>
+          }
+   finish {
+      // make it look wood-like
+      ambient 0.15
+      diffuse 0.6
+      // make it a little bit shiny
+      specular 0.3 roughness 0.01
+      phong 0.3 phong_size 60
+   }
+}
+
+
+// ttfn!

cdg0gv-2010/beispiele/triblend_off.p

@@ -0,0 +1,186 @@
+program triexample;
+uses geograph,hiddenl;
+
+type  funct2d = function(x,y: real) : real;   
+      funct2d2d = procedure(x,y: real; var gx,gy: real);
+      funct3d = function(p: vt3d) : real;   
+      funct3d3d = procedure(p: vt3d; var g: vt3d);
+      funct2d3d = procedure(u,v: real; var p: vt3d);
+      funct1d3d = procedure(u: real; var p: vt3d);
+      psurface_tangents = procedure(u,v: real; var su,sv: vt3d); 
+      implicit3d = procedure(p: vt3d; var fvalue: real; var gradf: vt3d);
+
+{var   f_gradf: implicit3d;}
+{--------------------------}
+{for triangulation:}
+const tnfmax=15000; tnemax=30000; tnpmax=15000; 
+      tfrnpmax= 500;  {max. number of points in a front}  
+      tnfrmax=20;     {max. number of additional fronts.}
+{--------------------------}
+var   i,ipro,ianf    : integer;
+      oriented_faces : boolean;
+      mue,dw,r_sph   : real;
+      p1	     : vt3d;
+
+{****************}
+{$i triang_proc.p}   
+{****************}
+
+{###################################################################}
+{****  auxiliary procs for generating OFFfiles used for GEOMVIEW ***}
+{*******************************************************************}
+
+{*****}
+var offdatei : text;
+{*****}
+   
+procedure open_offdatei;
+{ for OFF-file }
+ var  datei : string;
+ begin
+   writeln('name of OFF-File ? (... .off)');
+   readln(datei);
+   assign(offdatei,datei);
+   rewrite(offdatei);
+ end; { open_offdatei }
+{*************}
+
+procedure close_offdatei;
+ begin
+   Close(offdatei);
+ end;
+{************}
+
+procedure write_nangles_to_offfile;
+var i,k : integer;
+begin
+ open_offdatei;
+ writeln(offdatei, 'OFF');
+ writeln(offdatei,np,' ',nf,' ',0);
+ for i:= 1 to np do with p[i] do
+           writeln(offdatei,'   ',x:3:5,' ',y:3:5,' ',z:3:5);
+ for i:= 1 to nf do
+    with face[i] do
+         begin
+         write(offdatei,npf,'  ');
+         for k:= 1 to npf do write(offdatei,fp[k]-1,' ');
+         writeln(offdatei,' ');
+	 end;   
+ close_offdatei;
+end;  { write_nangles_to_offdatei }
+{#####################################################################}
+
+procedure f_gradf_blend3cy(p: vt3d; var f: real; var gradf: vt3d); 
+{Blending 3 cylinders by parabol. funct. splines}
+var c1,c2,c3,fk: real;
+begin
+  with p do 
+      begin
+      c1:= y*y+z*z-1;        c2:= x*x+z*z-1;        c3:= x*x+y*y-1;
+      fk:= 1+sqr(r_sph) - x*x - y*y - z*z;
+      f:= (1-mue)*c1*c2*c3-mue*fk*fk*fk;
+      gradf.x:= (1-mue)*2*x*c1*(c3 + c2) + 6*mue*x*fk*fk;
+      gradf.y:= (1-mue)*2*y*c2*(c3 + c1) + 6*mue*y*fk*fk;
+      gradf.z:= (1-mue)*2*z*c3*(c1 + c2) + 6*mue*z*fk*fk;
+      end;
+end;   
+{******}
+{*************************************************************}
+begin {main program}
+   graph_on(0);
+   repeat
+     writeln('****************************************************');
+     writeln('***     Triangulation of an implicit surface    *** ');   
+     writeln('***     (blending of 3 cylinders)               *** ');   
+     writeln('****************************************************'); 
+     writeln;
+     writeln('   and generation of a OFF-file for GEOMVIEW        ');  
+     writeln('****************************************************');  
+     writeln;
+
+{------------------------------------}      
+{for triangulation:}
+     cuttype:=0; {box}
+     {cuttype=1: cylinder, needs xcut,ycut,rcut_square (see proc. cut...)} 
+     {cuttype=2: sphere, needs xcut,ycut,zcut,rcut_square (   "         )} 
+     for i:= 1 to tnpmax do tpoint[i].full:= false;     
+     tnp:= 0;   tnf:= 0;    tnfr:= 0;
+       np:=0;    nf:= 0;      ne:= 0;   
+
+     writeln('Number of triangles ?');  readln(n_triang);
+     writeln;
+
+{------------------------------------------------------------------------}
+
+{Tringulation of a blend surface of 3 cylinders with bounding box}
+     f_gradf:= f_gradf_blend3cy;   mue:= 0.0003;  r_sph:= 3;  
+     put3d(1,1,-1, p1);       
+     dw:= pi2/30;   tstepl:= dw;  {radii= 1 !!}
+     xmin:= -3;  xmax:= 3;         
+     ymin:= -3;  ymax:= 3;
+     zmin:= -3;  zmax:= 3;
+     tnp:= 0;        tnf:= 0;    tnfr:= 0; 
+     start_triangulation(p1,f_gradf);
+
+{triangulation:}
+     nearpointtest:= true;
+     triangulation(nearpointtest);
+
+     writeln('tfrontnp: ',tfrontnp,' fullcount: ',fullcount,' tnf: ',
+             tnf,' n_triang: ',n_triang);
+	writeln;
+	writeln('tnp: ',tnp,'   tnf: ',tnf);  
+	writeln;
+
+{tri --> hiddenl.:}
+       for i:= 1 to tnp do p[np+i]:= tpoint[i].p;
+       for i:= 1 to tnf do 
+           with face[nf+i] do 
+                begin 
+                npf:= 3; 
+                with tface[i] do 
+                     begin fp[1]:= np+p1; fp[2]:= np+p2; fp[3]:= np+p3;  end; 
+                end; 
+	np:=np+tnp;  nf:= nf+tnf;
+       aux_polyhedron;   
+       writeln('np: ',np);  writeln('nf: ',nf); writeln('ne: ',ne); 
+
+{--------------------------------------------------------------------}
+     repeat
+     init_centralparallel_projection(2);
+
+{drawing : }
+     draw_area(250,250,120,120,25); 
+
+{hiddenline:}
+     oriented_faces:= false;
+     cp_lines_before_convex_faces(oriented_faces,true,false);
+  
+     draw_end;        writeln ;
+ 
+     writeln('Another projection?  (yes: 1, no: 0)');
+     readln(ipro);
+     until ipro=0;
+
+   writeln('Run again ?  (yes: 1, no: 0)');
+   readln(ianf);
+   until ianf=0;
+
+{generates offfile for GEOMVIEW:}
+   write_nangles_to_offfile;
+
+   graph_off;
+ end.
+
+
+
+
+
+
+
+
+
+
+
+
+

cdg0gv-2010/beispiele/triblend_pov.p

@@ -0,0 +1,196 @@
+program triexample;
+uses geograph,hiddenl;
+
+type  funct2d = function(x,y: real) : real;   
+      funct2d2d = procedure(x,y: real; var gx,gy: real);
+      funct3d = function(p: vt3d) : real;   
+      funct3d3d = procedure(p: vt3d; var g: vt3d);
+      funct2d3d = procedure(u,v: real; var p: vt3d);
+      funct1d3d = procedure(u: real; var p: vt3d);
+      psurface_tangents = procedure(u,v: real; var su,sv: vt3d); 
+      implicit3d = procedure(p: vt3d; var fvalue: real; var gradf: vt3d);
+
+{var   f_gradf: implicit3d;}
+{--------------------------}
+{for triangulation:}
+const tnfmax=15000; tnemax=30000; tnpmax=15000; 
+      tfrnpmax= 500;  {max. number of points in a front}  
+      tnfrmax=20;     {max. number of additional fronts.}
+{--------------------------}
+var   i,ipro,ianf    : integer;
+      oriented_faces : boolean;
+      mue,dw,r_sph   : real;
+      p1	     : vt3d;
+
+{****************}
+{$i triang_proc.p}   
+{****************}
+
+{****************************************************************}
+{****  auxiliary proc. for generating ppfiles used for POVRAY ***}
+{****************************************************************}
+   
+var ppfile : text;
+{*******}
+procedure open_ppfile;
+{ for ppfile }
+ var  datei : string;
+ begin
+   writeln('name of ppfile ? (... .pp)');
+   readln(datei);
+   assign(ppfile,datei);
+   rewrite(ppfile);
+ end; { open_ppfile }
+{*************}
+procedure close_ppdatei;
+ begin
+   Close(ppfile);
+ end;
+{************}
+procedure write_s3angles_to_ppfile(n1,n2: integer);
+var i : integer;
+begin
+ open_ppfile;
+ for i:= n1 to n2 do
+     begin
+     writeln(ppfile, 'smooth_triangle{');
+     with tface[i] do
+          begin
+          with tpoint[p1] do
+               begin
+               with p do
+                      write(ppfile,'<',x:3:5,', ',y:3:5,', ',z:3:5,'>,');
+               with nv do      
+                    writeln(ppfile,'<',x:3:5,', ',y:3:5,', ',z:3:5,'>,');
+               end;           
+          with tpoint[p2] do
+               begin
+               with p do
+                      write(ppfile,'<',x:3:5,', ',y:3:5,', ',z:3:5,'>,');
+               with nv do      
+                    writeln(ppfile,'<',x:3:5,', ',y:3:5,', ',z:3:5,'>,');
+               end;           
+          with tpoint[p3] do
+               begin
+               with p do
+                      write(ppfile,'<',x:3:5,', ',y:3:5,', ',z:3:5,'>,');
+               with nv do      
+                    writeln(ppfile,'<',x:3:5,', ',y:3:5,', ',z:3:5,'>');
+               end;           
+          writeln(ppfile,'}');
+          end; { with }
+     end;  { for }
+ close_ppdatei;
+end;  { write_3angles_to_ppfile }
+{***********}
+
+{************************************************************}
+procedure f_gradf_blend3cy(p: vt3d; var f: real; var gradf: vt3d); 
+{Blending 3 cylinders by parabol. funct. splines}
+var c1,c2,c3,fk: real;
+begin
+  with p do 
+      begin
+      c1:= y*y+z*z-1;        c2:= x*x+z*z-1;        c3:= x*x+y*y-1;
+      fk:= 1+sqr(r_sph) - x*x - y*y - z*z;
+      f:= (1-mue)*c1*c2*c3-mue*fk*fk*fk;
+      gradf.x:= (1-mue)*2*x*c1*(c3 + c2) + 6*mue*x*fk*fk;
+      gradf.y:= (1-mue)*2*y*c2*(c3 + c1) + 6*mue*y*fk*fk;
+      gradf.z:= (1-mue)*2*z*c3*(c1 + c2) + 6*mue*z*fk*fk;
+      end;
+end;   
+{******}
+{*************************************************************}
+begin {main program}
+   graph_on(0);
+   repeat
+     writeln('****************************************************');
+     writeln('***     Triangulation of implicit surfaces      *** ');   
+     writeln('****************************************************'); 
+     writeln;
+     writeln(' and generation of a PP-file for POVRAY (ray tracer) ');  
+     writeln('****************************************************');  
+     writeln;
+
+{------------------------------------}      
+{for triangulation:}
+     cuttype:=0; {box}
+     {cuttype=1: cylinder, needs xcut,ycut,rcut_square (see proc. cut...)} 
+     {cuttype=2: sphere, needs xcut,ycut,zcut,rcut_square (   "         )} 
+     for i:= 1 to tnpmax do tpoint[i].full:= false;     
+     tnp:= 0;   tnf:= 0;    tnfr:= 0;
+       np:=0;    nf:= 0;      ne:= 0;   
+
+     writeln('Number of triangles ?');  readln(n_triang);
+     writeln;
+
+{------------------------------------------------------------------------}
+
+{Tringulation of a blend surface of 3 cylinders with bounding box}
+     f_gradf:= f_gradf_blend3cy;   mue:= 0.0003;  r_sph:= 3;  
+     put3d(1,1,-1, p1);       
+     dw:= pi2/30;   tstepl:= dw;  {radii= 1 !!}
+     xmin:= -3;  xmax:= 3;         
+     ymin:= -3;  ymax:= 3;
+     zmin:= -3;  zmax:= 3;
+     tnp:= 0;        tnf:= 0;    tnfr:= 0; 
+     start_triangulation(p1,f_gradf);
+
+{triangulation:}
+     nearpointtest:= true;
+     triangulation(nearpointtest);
+
+     writeln('tfrontnp: ',tfrontnp,' fullcount: ',fullcount,' tnf: ',
+             tnf,' n_triang: ',n_triang);
+	writeln;
+	writeln('tnp: ',tnp,'   tnf: ',tnf);  
+	writeln;
+
+{tri --> hiddenl.:}
+       for i:= 1 to tnp do p[np+i]:= tpoint[i].p;
+       for i:= 1 to tnf do 
+           with face[nf+i] do 
+                begin 
+                npf:= 3; 
+                with tface[i] do 
+                     begin fp[1]:= np+p1; fp[2]:= np+p2; fp[3]:= np+p3;  end; 
+                end; 
+	np:=np+tnp;  nf:= nf+tnf;
+       aux_polyhedron;   
+       writeln('np: ',np);  writeln('nf: ',nf); writeln('ne: ',ne); 
+
+{--------------------------------------------------------------------}
+     repeat
+     init_centralparallel_projection(2);
+
+{drawing : }
+     draw_area(250,250,120,120,25); 
+
+{hiddenline:}
+     oriented_faces:= false;
+     cp_lines_before_convex_faces(oriented_faces,true,false);
+  
+     draw_end;        writeln ;
+ 
+     writeln('Another projection?  (yes: 1, no: 0)');
+     readln(ipro);
+     until ipro=0;
+
+   writeln('Run again ?  (yes: 1, no: 0)');
+   readln(ianf);
+   until ianf=0;
+
+{generates ppfile for POVRAY:}
+   write_s3angles_to_ppfile(1,tnf);
+
+   graph_off;
+ end.
+
+
+
+
+
+
+
+
+

cdg0gv-2010/beispiele/trisample.p

@@ -0,0 +1,342 @@
+program trisample;
+uses geograph,hiddenl;
+
+type  funct2d = function(x,y: real) : real;   
+      funct2d2d = procedure(x,y: real; var gx,gy: real);
+      funct3d = function(p: vt3d) : real;   
+      funct3d3d = procedure(p: vt3d; var g: vt3d);
+      funct2d3d = procedure(u,v: real; var p: vt3d);
+      funct1d3d = procedure(u: real; var p: vt3d);
+      psurface_tangents = procedure(u,v: real; var su,sv: vt3d); 
+      implicit3d = procedure(p: vt3d; var fvalue: real; var gradf: vt3d);
+
+{var   f_gradf: implicit3d;}
+{--------------------------}
+{for triangulation:}
+const tnfmax=15000; tnemax=30000; tnpmax=15000; 
+      tfrnpmax= 500;  {max. number of points in a front}  
+      tnfrmax=20;     {max. number of additional fronts.}
+{--------------------------}
+var   i,k,ipro,ianf,n1,n2,ik,ne0       		: integer;
+      u1,u2,v1,v2,du,dv,u,v			: real;
+      oriented_faces				: boolean;
+      x1,y1,r1,x2,y2,r2,x3,y3,r3,dw,r_tor,a_tor	: real;
+      p1					: vt3d;
+{****************}
+{$i triang_proc.p}   
+{****************}
+{************************************************************}
+
+procedure f_gradf_c1(p: vt3d; var f: real; var gradf: vt3d); 
+{1. cylinder}
+begin
+  with p do 
+  begin
+    f:= sqr(x-x1)+sqr(y-y1) -r1*r1;  
+    gradf.x:= 2*(x-x1);   
+    gradf.y:= 2*(y-y1);   
+    gradf.z:= 0;	 
+  end;
+end;   
+{******}
+
+procedure f_gradf_c2(p: vt3d; var f: real; var gradf: vt3d); 
+{2. cylinder}
+begin
+  with p do 
+  begin
+    f:= sqr(x-x2)+sqr(y-y2) -r2*r2;
+    gradf.x:= 2*(x-x2);   
+    gradf.y:= 2*(y-y2);   
+    gradf.z:= 0;	 
+  end;
+end;   
+{******}
+procedure f_gradf_sph4(p: vt3d; var f: real; var gradf: vt3d);  
+begin
+  with p do 
+       begin
+       f:= x*x*x*x+y*y*y*y+z*z*z*z - 16;
+       gradf.x:= 4*x*x*x;
+       gradf.y:= 4*y*y*y;
+       gradf.z:= 4*z*z*z;
+       end; 
+end;   
+{******}
+procedure f_gradf_tor(p: vt3d; var f: real; var gradf: vt3d); 
+{torus }
+var c1,c2: real;
+begin
+  c1:= sqr(r_tor);
+  with p do 
+      begin
+      c2:= x*x+y*y +z*z + c1 - sqr(a_tor);
+      f:= c2*c2 - 4*c1*(x*x+y*y);
+      gradf.x:= 4*x*c2 - 8*c1*x;
+      gradf.y:= 4*y*c2 - 8*c1*y;
+      gradf.z:= 4*z*c2;
+      end;
+end;   
+{******}
+{*************************************************************}
+begin {main program}
+   graph_on(0);
+   repeat
+     writeln('*********************************************');
+     writeln(' *** Triangulation of implicit surfaces  *** ');   
+     writeln('*********************************************'); 
+     writeln;
+     writeln('  example: 3 cylinders and a quartic sphere  ');  
+     writeln('*********************************************');  
+     writeln;
+{------------------------------------}      
+{for hiddenline:}
+      for i:= 1 to nemax do
+	  with edge[i] do begin color:=black; linew:=1; end;
+{------------------------------------}      
+{for triangulation:}
+     cuttype:=0; {box}
+     {cuttype=1: cylinder, needs xcut,ycut,rcut_square (see proc. cut...)} 
+     {cuttype=2: sphere, needs xcut,ycut,zcut,rcut_square (   "         )} 
+     for i:= 1 to tnpmax do tpoint[i].full:= false;     
+     tnp:= 0;   tnf:= 0;    tnfr:= 0;
+      np:=0;    nf:= 0;      ne:= 0;  ne0:= 0;
+
+     writeln('Number of triangles ?');  readln(n_triang);
+     writeln;
+
+{------------------------------------------------------------------------}
+
+{1. cylinder: with start hexagon}
+     f_gradf:= f_gradf_c1;
+     x1:= -2; y1:= -5;	 r1:= 2; 
+     put3d(0,-3,1, p1);   
+     tstepl:= 0.4; 
+     xmin:= -10;  xmax:= 10;         
+     ymin:= -10;  ymax:= 10;
+     zmin:= -2;   zmax:= 5;
+     tnp:= 0;        tnf:= 0;    tnfr:= 0; 
+     start_triangulation(p1,f_gradf); 
+
+{triangulation: }
+     nearpointtest:= true;
+     triangulation(nearpointtest);
+
+     writeln('tfrontnp: ',tfrontnp,' fullcount: ',fullcount,' tnf: ',
+             tnf,' n_triang: ',n_triang);
+     writeln;
+     writeln('tnp: ',tnp,'   tnf: ',tnf);  
+     writeln;
+
+{tri --> hiddenl.:}
+     for i:= 1 to tnp do p[np+i]:= tpoint[i].p;
+     for i:= 1 to tnf do 
+         with face[nf+i] do 
+              begin 
+              npf:= 3; 
+              with tface[i] do 
+	      begin fp[1]:= np+p1; fp[2]:= np+p2; fp[3]:= np+p3;  end;
+	      end; 
+     np:=np+tnp;  nf:= nf+tnf;
+     aux_polyhedron;   
+     writeln('np: ',np);  writeln('nf: ',nf); writeln('ne: ',ne); 
+     for i:= ne0+1 to ne do
+	  with edge[i] do begin color:=black; linew:=1; end;
+      ne0:= ne;
+{--------------------------------}
+
+{2. cylinder: with start POLYGON (circle)}
+     for i:= 1 to tnpmax do tpoint[i].full:= false;     
+     tnp:= 0;        tnf:= 0;    tnfr:= 0;
+     x2:= -5;  y2:=0;  r2:= 2;
+     f_gradf:= f_gradf_c2;   
+     xmin:= -10;  xmax:= 10;         
+     ymin:= -10;  ymax:= 10;
+     zmin:= -10;  zmax:= 10;
+     tfrontnp:= 30; tnp:= tfrontnp;  dw:= pi2/tfrontnp;   tstepl:= r2*pi2/tfrontnp;
+	for i:= 1 to tfrontnp do   {front polygon: circle on the top} 
+         begin
+         with tpoint[i] do 
+              begin
+              put3d(x2+r2*cos((i-1)*dw),y2+r2*sin((i-1)*dw),5, p1); 
+              surface_point_normal_tangentvts(p1,f_gradf, p,nv,tv1,tv2);
+              full:= false; achange:= true;
+              end;  { with }
+         tfrontpt[i]:= i;
+         end;  { for }  
+     dw:= -dw;    
+	for i:= 1 to tfrontnp do    {bounding polygon: circle on the bottom}
+         begin
+         with tpoint[tnp+i] do 
+              begin
+              put3d(x2+r2*cos((i-1)*dw),y2+r2*sin((i-1)*dw),-2, p1); 
+              surface_point_normal_tangentvts(p1,f_gradf, p,nv,tv1,tv2);
+              full:= false; achange:= true;
+              end;  { with }
+         tfr[1,i]:= tnp+i ;
+         end;  { for }
+     tfr[1,0]:= tfrontnp;   tnp:= tnp+tnp;    
+     tnfr:= 1;
+
+{triangulation:}
+     nearpointtest:= true;
+     triangulation(nearpointtest);
+
+     writeln('tfrontnp: ',tfrontnp,' fullcount: ',fullcount,' tnf: ',
+             tnf,' n_triang: ',n_triang);
+	writeln;
+	writeln('tnp: ',tnp,'   tnf: ',tnf);  
+	writeln;
+
+{tri --> hiddenl.:}
+       for i:= 1 to tnp do p[np+i]:= tpoint[i].p;
+       for i:= 1 to tnf do 
+           with face[nf+i] do 
+                begin 
+                npf:= 3; 
+                with tface[i] do 
+                     begin fp[1]:= np+p1; fp[2]:= np+p2; fp[3]:= np+p3;  end; 
+                end; 
+	np:=np+tnp;  nf:= nf+tnf;
+       aux_polyhedron;   
+       writeln('np: ',np);  writeln('nf: ',nf); writeln('ne: ',ne); 
+     for i:= ne0+1 to ne do
+	  with edge[i] do begin color:=red; linew:=3; end;
+      ne0:= ne;
+{-------------------------------------------------------}
+
+{3. quartic sphere x^4+y^4+z^4-16=0: with start hexagon}
+     for i:= 1 to tnpmax do tpoint[i].full:= false;     
+     tnp:= 0;        tnf:= 0;    tnfr:= 0;
+     f_gradf:= f_gradf_sph4;   
+     put3d(0,0,1, p1);
+     tstepl:= 0.3;                  
+     xmin:= -3;  xmax:= 3;         
+     ymin:= -3;  ymax:= 3;
+     zmin:= -3;  zmax:= 3;
+     tnp:= 0;        tnf:= 0;    tnfr:= 0; 
+     start_triangulation(p1,f_gradf); 
+
+{triangulation: }
+     nearpointtest:= true;
+     triangulation(nearpointtest);
+
+     writeln('tfrontnp: ',tfrontnp,' fullcount: ',fullcount,' tnf: ',
+             tnf,' n_triang: ',n_triang);
+     writeln;
+     writeln('tnp: ',tnp,'   tnf: ',tnf);  
+     writeln;
+
+{tri --> hiddenl.:}
+     for i:= 1 to tnp do p[np+i]:= tpoint[i].p;
+     for i:= 1 to tnf do 
+         with face[nf+i] do 
+              begin 
+              npf:= 3; 
+              with tface[i] do 
+	      begin fp[1]:= np+p1; fp[2]:= np+p2; fp[3]:= np+p3;  end;
+	      end; 
+     np:=np+tnp;  nf:= nf+tnf;
+     aux_polyhedron;   
+     writeln('np: ',np);  writeln('nf: ',nf); writeln('ne: ',ne); 
+     for i:= ne0+1 to ne do
+	  with edge[i] do begin color:=blue; linew:=1; end;
+      ne0:= ne;
+
+{--------------------------------}
+      
+{4.torus: with start hexagon}
+     for i:= 1 to tnpmax do tpoint[i].full:= false;     
+     tnp:= 0;        tnf:= 0;    tnfr:= 0;
+     f_gradf:= f_gradf_tor; r_tor:= 8.5; a_tor:= 1;   
+{     cuttype:= 1; xcut:= 0;  ycut:= 0; rcut_square:= sqr(r_tor);}
+     put3d(0,6,0, p1);
+     tstepl:= 0.45;                  
+     xmin:= -10;  xmax:= 10;         
+     ymin:= -10;  ymax:= 10;
+     zmin:= -3;  zmax:= 3;
+     start_triangulation(p1,f_gradf); 
+
+{triangulation: }
+     nearpointtest:= true;
+     triangulation(nearpointtest);
+
+     writeln('tfrontnp: ',tfrontnp,' fullcount: ',fullcount,' tnf: ',
+             tnf,' n_triang: ',n_triang);
+     writeln;
+     writeln('tnp: ',tnp,'   tnf: ',tnf);  
+     writeln;
+
+{tri --> hiddenl.:}
+     for i:= 1 to tnp do p[np+i]:= tpoint[i].p;
+     for i:= 1 to tnf do 
+         with face[nf+i] do 
+              begin 
+              npf:= 3; 
+              with tface[i] do 
+	      begin fp[1]:= np+p1; fp[2]:= np+p2; fp[3]:= np+p3;  end;
+	      end; 
+     np:=np+tnp;  nf:= nf+tnf;
+     aux_polyhedron;   
+     writeln('np: ',np);  writeln('nf: ',nf); writeln('ne: ',ne); 
+     for i:= ne0+1 to ne do
+	  with edge[i] do begin color:=green; linew:=1; end;
+      ne0:= ne;
+
+{--------------------------------}
+      
+{5. cylinder: with QUADRANGLES}
+     x3:= -2; y3:= 5; r3:= 2;
+     n1:= 40;    n2:= 10;
+     u1:= 0;     u2:= pi2;
+     v1:=-2;     v2:= 5;
+     du:= (u2-u1)/n1;     dv:= (v2-v1)/(n2-1);
+     v:= v1;
+     for k:= 1 to n2 do
+         begin
+         u:= u1;
+         for i:= 1 to n1 do
+             begin
+             ik:=i + (k-1)*n1;
+             put3d(x3+r3*cos(u),y3+r3*sin(u),v,  p[np+ik]);
+             u:= u+du;
+             end;
+         v:= v+dv;
+         end;
+     aux_cylinder(n1,n2,np,ne,nf);
+     for i:= ne0+1 to ne do
+	  with edge[i] do begin color:=cyan; linew:=1; end;
+      ne0:= ne;
+	
+{--------------------------------------------------------------------}
+     repeat
+     init_centralparallel_projection(2);
+
+{ for drawing : }
+     draw_area(250,250,120,120,10); 
+
+{hiddenline:}
+     oriented_faces:= false;   is_permitted:= true;    newstyles:= true;
+     cp_lines_before_convex_faces(oriented_faces,is_permitted,newstyles);
+  
+     draw_end;        writeln ;
+ 
+     writeln('Another projection?  (yes: 1, no: 0)');
+     readln(ipro);
+     until ipro=0;
+
+   writeln('Run again ?  (yes: 1, no: 0)');
+   readln(ianf);
+   until ianf=0;
+
+   graph_off;
+ end.
+
+
+
+
+
+
+
+
+

cdg0gv-2010/include/geoconst.pas

@@ -0,0 +1,9 @@
+{geoconst.pas:}
+   array_size = 10000; {  < array_size in driver/g_const.h !!!!}
+   pi= 3.14159265358;       pi2= 6.2831853;       pih= 1.5707963;
+   eps1=0.1;      eps2=0.01;       eps3=0.001;     eps4=0.0001;
+   eps5=0.00001;  eps6=0.000001;   eps7=0.0000001; eps8=0.00000001;
+   default=-1; black=0; blue=1; green=2; cyan=3; red=4; magenta=5; brown=6; 
+   lightgray=7; darkgray=8; lightblue=9; lightgreen=10; lightcyan=11; 
+   lightred=12; lightmagenta=13; yellow=14; white=15;
+

cdg0gv-2010/include/geoproc.pas

@@ -0,0 +1,232 @@
+{##############################}
+{ geoproc.pas for FREE-PASCAL }
+{##############################}
+
+const
+  GeoPLDcommand    : Boolean = false; 
+  GeoPLDinteraktiv : Boolean = false;
+  GeoBild          : Boolean = true;
+  GeoMono          : Boolean = false;
+
+var
+  GeoPLDFilename  : string;
+  GeoEpsFilename  : string;                             
+
+{*****************************************************}
+{*** A: For Generation of a PLD-file ***}
+{*****************************************************}
+{ i include/geoinit.pas}
+{$i pld.pas}
+{$i postscr.pas}
+
+{***********************************************************************}
+   {***  B: for the generation of an EPS-file of the drawing and ....***}
+{***********************************************************************}
+
+procedure graph_on(mode : integer);
+begin
+  null2d.x:= 0;  null2d.y:= 0;
+  null3d.x:= 0;  null3d.y:= 0;   null3d.z:= 0;
+(*   if GeoBild then
+   begin
+      GeoEpsFilename:= 'geodummy.eps';
+      ps_graph_on(GeoEpsFilename, true);
+   end; *)
+  GeoPLDinteraktiv := (mode <> 0);
+  plot_it := GeoPLDinteraktiv or GeoPLDcommand;
+  if plot_it then  plot_init;
+end; { graph_on }
+{*************}
+
+procedure draw_area(width,height,x0,y0,sfac : real);
+{.... and opens geodummy.eps}
+begin
+  if plot_it then  plot_area(width, height);
+  scalefactor:= sfac;
+  origin2d.x:=x0; origin2d.y:=y0;
+  if scalefactor<>1 then writeln('scalefactor is:',scalefactor:3:2);
+  if GeoBild then
+     begin
+      GeoEpsFilename:= 'geodummy.eps';
+      ps_graph_on(GeoEpsFilename, true);
+      ps_draw_area(width,height,x0,y0, true);                   
+     end;	
+ end; { draw_area}
+{*************}
+
+procedure draw_end;
+{closes the eps-file and displays it}
+var s :  longint;
+begin
+ if plot_it then plot_end;
+ if GeoBild then
+    begin
+       ps_draw_end;
+       writeln('###: show drawing (with gv):');
+       writeln('###: (continue after quitting gv)');
+       s:=fpsystem('gv geodummy.eps');
+{       writeln('shell beendet mit: ',s);}
+    end;
+end; { draw_end }
+{*************}
+
+procedure graph_off;
+begin
+  if GeoBild then
+     begin
+	ps_graph_off;
+     end;
+end; {graph_off}
+{**************}
+
+procedure new_color(color: integer);
+begin
+  if not GeoMono then 
+  begin
+    if plot_it then plot_new_color(color);
+    if color = default then color := black;
+    if GeoBild then
+       begin
+	  ps_new_color(color);                         
+       end;
+  end;
+end;  { new_color }
+{**************}
+
+procedure new_linewidth(factor: real);
+const
+   default = 3.0;
+begin
+   if plot_it then plot_new_linewidth(factor);
+   if GeoBild then
+      begin
+	 ps_new_linewidth(factor);          
+      end;
+end;
+
+procedure linec2d(x1,y1,x2,y2 : real; style : integer);
+{draws the line (x1,y1)(x2,y2)}
+begin
+  x1 := x1*scalefactor;
+  x2 := x2*scalefactor;
+  y1 := y1*scalefactor;
+  y2 := y2*scalefactor;
+  if style > 3 then style := 0;
+  if plot_it then plot_linec2d(x1,y1, x2,y2,style);
+  if GeoBild then
+     begin
+	ps_linec2d(x1,y1,x2,y2, style);                  
+     end;
+end; {linec2d}
+{*************}
+
+procedure line2d(p1,p2 : vt2d; style : integer);
+begin
+  linec2d(p1.x,p1.y,p2.x,p2.y,style);
+end; {line2d}
+{*************}
+
+procedure pointc2d(x,y : real; style : integer);
+{draws a point}
+begin
+  x:= x*scalefactor;
+  y:= y*scalefactor;
+  if plot_it  then plot_pointc2d(x,y,style);
+  if GeoBild then
+     begin
+	ps_pointc2d(x,y,style);
+     end;
+end; { pointc2d }
+{*************}
+
+procedure point2d(p : vt2d; style: integer);
+begin
+  pointc2d(p.x,p.y,style);
+end;  {point2d}
+{**************}
+
+procedure curve2d(var p: vts2d; n1,n2,style : integer);
+{draws the polyline (x(n1),y(n1))... (x(n2),y(n2))}
+var
+  i  : integer;
+  pp : vts2d;
+begin
+    for i:= n1 to n2 do 
+    begin
+      pp[i].x := p[i].x * scalefactor;
+      pp[i].y := p[i].y * scalefactor;
+    end;
+    if style > 3 then style := 0;
+    if plot_it then plot_curve2d(pp,n1,n2,style);
+    if GeoBild then
+       begin
+	  ps_curve2d(pp,n1,n2,style);
+       end;
+end; {curve2d}
+{*************}
+
+procedure curve2d_vis(var p: vts2d; n1,n2,style: integer; visible: b_array );
+{neighbored "visible" points are connected.
+ if style=10: "invisible" lines are dashed. }
+ var i,i1 : integer;
+     vis : boolean;
+ begin
+
+   i1 := n1;
+   i := i1;
+   vis := visible[i1];
+   
+   repeat
+     i := i+1;
+     if (i>n2) or (vis <> visible[i]) then 
+     begin
+        if vis then curve2d(p,i1,i-1,style)
+        else if style=10 then curve2d(p,i1,i-1,1);
+        i1 := i;
+        vis := visible[i];
+     end;  
+   until i>n2;
+
+ end; { curve2d_vis }
+{*************}
+
+procedure arrowc2d(x1,y1,x2,y2 : real; style : integer);
+{draws an arrow}
+ var x21,y21,x3,y3,x4,y4,d,sl,sb,sld,sbd : real;
+ begin
+   x21:= x2-x1;   y21:= y2-y1;
+   d:= sqrt(x21*x21+y21*y21);
+   sl:= 3/scalefactor;  {Laenge der Spitze}
+   sb:= 1/scalefactor; {Breite ..}
+   if d>=sl then   {Pfeilspitze}
+              begin
+                sld:= sl/d;  sbd:= sb/d;
+                x3:= x2-sld*x21-sbd*y21 ;   y3:= y2-sld*y21+sbd*x21 ;
+                x4:= x2-sld*x21+sbd*y21 ;   y4:= y2-sld*y21-sbd*x21 ;
+                linec2d(x2,y2,x3,y3,0)  ;   linec2d(x2,y2,x4,y4,0)  ;
+              end;
+   linec2d(x1,y1,x2,y2,style);
+ end; { arrowc2d }
+{**************}
+
+procedure arrow2d(p1,p2 : vt2d; style : integer);
+ begin
+   arrowc2d(p1.x,p1.y,p2.x,p2.y,style);
+ end; { arrow2d }
+{**********************************************************************}
+
+procedure read_integer_file(file_name: string; n_dat: integer;
+                                           var int_var: i_array);
+var  text_var : text;
+     i   : integer;
+begin
+   assign(text_var,file_name);
+   reset(text_var);
+   for i:= 1 to n_dat do read(text_var, int_var[i]);
+   close(text_var);  
+end;  { read_integer_file }
+{***********}
+
+
+
+

cdg0gv-2010/include/geotype.pas

@@ -0,0 +1,10 @@
+{geotype.pas:}
+
+       r_array = array [0..array_size] of real;
+       i_array = array [0..array_size] of integer;
+       b_array = array [0..array_size] of boolean;
+       vt2d    = record  x: real;  y: real; end;
+       vt3d    = record  x: real;  y: real;  z: real; end;
+       vts2d   = array[0..array_size] of vt2d;
+       vts3d   = array[0..array_size] of vt3d;
+       matrix3d= array[1..3,1..3] of real;

cdg0gv-2010/include/geovar.pas

@@ -0,0 +1,19 @@
+ {geovar.pas:}
+
+       null2d : vt2d;  null3d : vt3d;            {Nullvektoren}
+       scalefactor: real;                        
+
+    {**fuer area_2d and curve2d:}
+       origin2d  : vt2d;
+
+    {**fuer Parallel- und Zentral-Projektion:}
+       u_angle,v_angle,                     {Projektionswinkel}
+       rad_u,rad_v,                         {rad(u), rad(v)}
+       sin_u,cos_u,sin_v,cos_v : real;      {sin-,cos- Werte von u, v}
+       e1vt,e2vt,n0vt : vt3d;               {Basis-Vektoren}
+                                            {Normalen-Vektor der Bildebene}
+    {**fuer Zentral-Projektion:}
+       mainpt,                              {Hauptpunkt}
+       centre : vt3d;                       {Zentrum}
+       distance : real;                     {Distanz Hauptpunkt-Zentrum}
+

cdg0gv-2010/include/head_ag.pas

@@ -0,0 +1,75 @@
+{**********************************************}
+{********* H E A D _ A G  *********************}
+{*** Prozedur-Koepfe der Datei proc_ag.pas  ***}
+{**********************************************}
+
+function sign(a: real) : integer;
+function max(a,b : real) : real;
+function min(a,b : real) : real;
+procedure change1d(var a,b: real);
+procedure change2d(var v1,v2: vt2d);
+procedure change3d(var v1,v2: vt3d);
+procedure put2d(x,y : real; var v: vt2d);
+procedure put3d(x,y,z : real; var v: vt3d);
+procedure get3d(v : vt3d; var x,y,z: real);
+procedure scale2d(r : real; v: vt2d; var vs: vt2d);
+procedure scale3d(r : real; v: vt3d; var vs: vt3d);
+procedure scaleco2d(r1,r2 : real; v: vt2d; var vs: vt2d);
+procedure scaleco3d(r1,r2,r3 : real; v: vt3d; var vs: vt3d);
+procedure sum2d(v1,v2 : vt2d; var vs : vt2d);
+procedure sum3d(v1,v2 : vt3d; var vs : vt3d);
+procedure diff2d(v1,v2 : vt2d; var vs : vt2d);
+procedure diff3d(v1,v2 : vt3d; var vs : vt3d);
+procedure lcomb2vt2d(r1: real; v1: vt2d; r2: real; v2: vt2d; var vlc : vt2d);
+procedure lcomb2vt3d(r1: real; v1: vt3d; r2: real; v2: vt3d; var vlc : vt3d);
+procedure lcomb3vt2d(r1: real; v1: vt2d; r2: real; v2: vt2d;
+                                         r3: real; v3: vt2d; var vlc : vt2d);
+procedure lcomb3vt3d(r1: real; v1: vt3d; r2: real; v2: vt3d;
+                                          r3: real; v3: vt3d; var vlc : vt3d);
+procedure lcomb4vt3d(r1: real; v1: vt3d; r2: real; v2: vt3d;
+                     r3: real; v3: vt3d; r4: real; v4: vt3d; var vlc : vt3d);
+function abs2d(v : vt2d) : real;
+function abs3d(v : vt3d) : real;
+function length2d(v : vt2d) : real;
+function length3d(v : vt3d) : real;
+function scalarp2d(p1,p2 : vt2d) : real;
+function scalarp3d(p1,p2 : vt3d) : real;
+procedure normalize2d(var p: vt2d);
+procedure normalize3d(var p: vt3d);
+function distance2d(p,q : vt2d): real;
+function distance3d(p,q : vt3d): real;
+function distance2d_square(p,q : vt2d) : real;
+function distance3d_square(p,q : vt3d) : real;
+procedure vectorp(v1,v2 : vt3d; var vp : vt3d);
+function determ3d(v1,v2,v3: vt3d) : real;
+procedure rotor2d(cos_rota,sin_rota : real; p : vt2d; var pr: vt2d);
+procedure rotp02d(cos_rota,sin_rota : real; p0,p : vt2d; var pr: vt2d);
+procedure rotorz(cos_rota,sin_rota : real; p : vt3d; var pr: vt3d);
+procedure rotp0z(cos_rota,sin_rota : real; p0,p: vt3d; var pr: vt3d);
+procedure rotp0x(cos_rota,sin_rota : real; p0,p: vt3d; var pr: vt3d);
+procedure rotp0y(cos_rota,sin_rota : real; p0,p: vt3d; var pr: vt3d);
+function polar_angle(x,y : real) : real;
+procedure equation_degree1(a,b :real; var x : real ; var ns : integer);
+procedure equation_degree2(a,b,c:real; var x1,x2:real ;var ns:integer);
+procedure is_line_line(a1,b1,c1, a2,b2,c2 : real; var xs,ys : real;
+                                                  var nis : integer);
+procedure is_unitcircle_line(a,b,c: real; var x1,y1,x2,y2 : real;
+                                                  var nis : integer);
+procedure is_circle_line(xm,ym,r, a,b,c: real; var x1,y1,x2,y2 : real;
+                                               var nis : integer);
+procedure is_circle_circle(xm1,ym1,r1,xm2,ym2,r2: real;
+                                 var x1,y1,x2,y2: real; var nis: integer);
+function pt_before_plane(p,nv: vt3d; d: real) : boolean;
+procedure plane_equ(p1,p2,p3 : vt3d; var nv : vt3d; var d: real;
+                                                      var error : boolean);
+procedure is_line_plane(p,rv,nv: vt3d; d : real; var pis : vt3d;
+                                                 var nis : integer);
+procedure is_3_planes(nv1: vt3d; d1: real; nv2: vt3d; d2: real;
+                      nv3: vt3d; d3: real;
+                      var pis : vt3d; var error: boolean);
+procedure is_plane_plane(nv1: vt3d; d1: real; nv2: vt3d; d2: real;
+                             var p,rv : vt3d; var error: boolean);
+procedure ptco_plane3d(p0,v1,v2,p: vt3d; var xi,eta: real; var error: boolean);
+procedure newcoordinates3d(p,b0,b1,b2,b3: vt3d; var pnew: vt3d);
+
+

cdg0gv-2010/include/head_geo.pas

@@ -0,0 +1,21 @@
+{***************************}
+{***  H E A D _ G E O    ***}
+{***************************}
+
+(*procedure GeoInit;*)
+procedure graph_on(mode : integer);
+procedure draw_area(width,height,x0,y0,sfac : real);
+procedure draw_end;
+procedure graph_off;
+procedure new_color(color: integer);
+procedure new_linewidth(factor :real);
+procedure linec2d(x1,y1,x2,y2 : real; style : integer);
+procedure line2d(p1,p2 : vt2d; style : integer);
+procedure pointc2d(x,y : real; style : integer);
+procedure point2d(p : vt2d; style: integer);
+procedure curve2d(var p: vts2d; n1,n2,style : integer);
+procedure curve2d_vis(var p: vts2d; n1,n2,style: integer; visible: b_array );
+procedure arrowc2d(x1,y1,x2,y2 : real; style : integer);
+procedure arrow2d(p1,p2 : vt2d; style : integer);
+procedure read_integer_file(file_name: string; n_dat: integer;
+                                           var int_var: i_array);

cdg0gv-2010/include/head_pp.pas

@@ -0,0 +1,17 @@
+{***********************}
+{***  H E A D _ P P  ***}
+{***********************}
+
+procedure init_parallel_projection;
+procedure pp_vt3d_vt2d(p : vt3d; var pp: vt2d);
+procedure pp_point(p: vt3d; style: integer);
+procedure pp_line(p1,p2 : vt3d ; style : integer);
+procedure pp_arrow(p1,p2 : vt3d; style : integer);
+procedure pp_axes(al : real);
+procedure pp_vts3d_vts2d(var p: vts3d; n1,n2 : integer; var pp : vts2d);
+procedure pp_curve(var p: vts3d; n1,n2,style : integer);
+procedure pp_curve_vis(var p : vts3d; n1,n2,style : integer; visible: b_array );
+procedure pp_line_before_plane(p1,p2,nv: vt3d; d : real; side,style: integer);
+procedure pp_curve_before_plane(var p: vts3d ; n1,n2: integer; nv: vt3d; d: real;
+                                                       side,style : integer);
+

cdg0gv-2010/include/head_zp.pas

@@ -0,0 +1,48 @@
+{*******************************}
+{*****  H E A D _ Z P      *****}
+{*****  Zentralprojektion  *****}
+{*****  OHNE Clipping      *****}
+{*******************************}
+
+procedure init_central_projection;
+{**************}
+
+procedure init_centralparallel_projection(ind : integer);
+{**************}
+
+procedure transf_to_e1e2n0_base(p : vt3d; var pm : vt3d);
+{***************}
+
+procedure cp_vt3d_vt2d(p: vt3d; var pp : vt2d);
+{**************}
+
+procedure cp_point(p: vt3d; style: integer);
+{*************}
+
+procedure cp_line(p1,p2 : vt3d ; style : integer);
+{*************}
+
+procedure cp_arrow(p1,p2 : vt3d; style : integer);
+{*************}
+
+procedure cp_axes(al : real);
+{*************}
+
+procedure cp_vts3d_vts2d(var p: vts3d; n1,n2 : integer; var pp : vts2d);
+{*************}
+
+procedure cp_curve(var p: vts3d; n1,n2,style : integer);
+{*************}
+
+
+ procedure cp_curve_vis(var p : vts3d; n1,n2,style : integer; visible: b_array );
+{*************}
+
+procedure cp_line_before_plane(p1,p2,nv: vt3d; d : real; side,style: integer);
+{**************}
+
+procedure cp_curve_before_plane(var p: vts3d ; n1,n2: integer; nv: vt3d; d: real;
+                                                       side,style : integer);
+{*************}
+
+

cdg0gv-2010/include/head_zpo.pas

@@ -0,0 +1,27 @@
+{*****************}
+{***  HEAD_ZPO ***}
+{*****************}
+
+procedure aux_polyhedron;
+procedure cp_vts3d_vts2d_spez(var p: vts3d; n1,n2: integer;
+                              var pp: vts2d; var pdist: r_array);
+procedure aux_quadrangle(n1,n2,np0,ne0,nf0: integer);
+procedure aux_quadrangle_triang(n1,n2: integer;
+                                show_triangles: boolean);
+procedure aux_cylinder(n1,n2,np0,ne0,nf0: integer);
+procedure aux_torus(n1,n2,np0,ne0,nf0: integer);
+procedure is_line_convex_polygon(p1,p2 : vt2d; p_pol : vts2d_pol; np : integer;
+                                           var t1,t2 : real; var ind : integer);
+procedure intmint(a,b,c,d: real; var e1,f1,e2,f2: real; var ind: integer);
+procedure cp_lines_before_convex_faces(oriented_faces,is_permitted,newstyles : boolean);
+procedure is_interv_interv(var a,b,c,d,aa,bb : real; var inters: boolean);
+procedure box3d_of_pts(var p : vts3d_pol; np: integer; var box : box3d_dat);
+function is_two_boxes3d(var box1,box2 : box3d_dat) : boolean;
+procedure is_line_conv_pol_in_plane3d(var pl,rl: vt3d; var pp : vts3d_pol;
+                                      npp : integer;
+                                      var t1,t2 : real; var inters : boolean);
+procedure is_n1gon_n2gon3d(var pp1,pp2: vts3d_pol; np1,np2: integer;
+                           var  ps1,ps2 : vt3d; var intersection : boolean);
+procedure boxes_of_faces;
+procedure is_face_face(i,k: integer;  var ps1,ps2 : vt3d;
+                                      var intersection: boolean);

cdg0gv-2010/include/pld.pas

@@ -0,0 +1,99 @@
+{pld.pas}
+
+{******************************************************}
+{*** A: Procedures for the generation of a PLD-file ***}
+{******************************************************}
+
+const
+      pt_per_mm : real = 10.0;
+      pld_ver   = 'Ver1.0';
+
+var
+      plotdatei : text;
+      plot_it   : boolean;
+
+procedure plot_init;
+ var i : integer;
+ begin
+   writeln('*************************************************************');
+   writeln('** A PLD-file (point-line-description) will be generated ****');
+   writeln('*************************************************************');
+   writeln;
+   writeln;
+ end; { plot_init }
+
+{*************}
+
+procedure plot_area(width,height: real);
+{ for PLD-file }
+ var  datei : string;
+ begin
+   if GeoPLDinteraktiv then begin
+     writeln('name of PLD-file ? (... .pld)');
+     readln(datei);
+   end;
+   if GeoPLDcommand then datei:=GeoPLDFilename;
+   assign(plotdatei,datei);
+   rewrite(plotdatei);
+   writeln(plotdatei,pld_ver);
+   Write(plotdatei,width:10:3,' ',height:10:3,' ');
+ end; { plot_area }
+{*************}
+
+procedure plot_end;
+ begin
+   Write(plotdatei,'*');
+   Close(plotdatei);
+ end;
+{*************}
+
+procedure plot_linec2d(x1,y1,x2,y2 : real; style : integer);
+{draws the line (x1,y1)(x2,y2)}
+ var px1,py1,px2,py2 : longint; 
+ begin
+   px1:=round((x1+origin2d.x)*pt_per_mm); 
+   py1:=round((y1+origin2d.y)*pt_per_mm);
+   px2:=round((x2+origin2d.x)*pt_per_mm); 
+   py2:=round((y2+origin2d.y)*pt_per_mm);
+   if style > 3 then style := 0;
+   Write(plotdatei,'L ',px1,' ',py1,' ',px2,' ',py2,' ',style,' ');
+ end; { plot_linec2d }
+{*************}
+
+procedure plot_curve2d(var p: vts2d; n1,n2,style : integer);
+{draws the polyline (x(n1),y(n1))... (x(n2),y(n2))}
+var
+  i : integer;
+  px1,py1 : longint;
+begin
+  if style > 3 then style := 0;
+  Write(plotdatei,'K ',n2-n1+1,' ',style,' ');
+  for i := n1 to n2 do
+  begin
+    px1:=round((p[i].x+origin2d.x)*pt_per_mm);    
+    py1:=round((p[i].y+origin2d.y)*pt_per_mm);
+    Write(plotdatei,px1,' ',py1,' ');
+  end;
+end; {plot_curve2d}
+{*************}
+
+procedure plot_pointc2d(x,y : real; style : integer);
+{draws a points}
+ var ix,iy : longint;
+ begin
+   ix:=round((x+origin2d.x)*pt_per_mm);
+   iy:=round((y+origin2d.y)*pt_per_mm);
+   Write(plotdatei, 'P ',ix,' ',iy,' ',style,' ')
+ end; { plot_pointc2d}
+{*************}
+
+procedure plot_new_color(color:integer);
+begin
+   Writeln(plotdatei, 'C',color)
+ end; { plot_new_color}
+{*************}
+
+procedure plot_new_linewidth(factor : real);
+begin
+   Writeln(plotdatei, 'W', round(factor*10));
+end;

cdg0gv-2010/include/postscr.pas

@@ -0,0 +1,228 @@
+{postscr.pas:}
+
+const
+   ps_style : Integer = 99;
+
+var
+   psfile : text;
+
+
+procedure ps_set_style(style :integer);
+begin
+  if (style <> ps_style) then
+  begin
+    Case style of
+          0  : writeln(psfile,'[] 0 setdash');
+          1  : writeln(psfile,'[2 1.5] 0 setdash');
+          3  : writeln(psfile,'[.2 1.5] 0 setdash');
+          2  : writeln(psfile,'[3 1 .2 1] 0 setdash');
+    end;
+    ps_style := style;
+  end;
+end;
+
+
+procedure ps_graph_on(fname:string;eps:boolean);
+begin
+  assign(psfile,fname);  
+  rewrite(psfile);
+  if not eps then
+  begin
+    writeln(psfile, '%!PS-Adobe-2.0 PSF-2.0');
+    writeln(psfile, '%%Creator: pldv');
+    writeln(psfile, '/Rahmen{newpath');
+    writeln(psfile, 'wx neg wy neg moveto ');
+    writeln(psfile, '0 hoehe     rlineto  ');
+    writeln(psfile, 'breite 0    rlineto  ');
+    writeln(psfile, '0 hoehe neg rlineto  ');
+    writeln(psfile, 'closepath            ');
+    writeln(psfile, '/lw currentlinewidth def ');
+    writeln(psfile, '0 setlinewidth           ');
+    writeln(psfile, 'stroke                   ');
+    writeln(psfile, 'lw setlinewidth          ');
+    writeln(psfile, 'wx neg wy neg 5 sub moveto ');
+    writeln(psfile, 'titel show             ');
+(*    writeln(psfile, '/seite seite 1 add def ');
+    writeln(psfile, 'wx 20 sub wy 5 add neg moveto ');
+    writeln(psfile, '(Seite ) show  ');
+    writeln(psfile, 'seite str cvs  ');*)
+    writeln(psfile, 'str show       ');
+    writeln(psfile, '}def           ');
+    writeln(psfile, '/titel (',fname ,') def');
+  end else {EPS}
+  begin
+    writeln(psfile, '%!PS-Adobe-2.0 EPSF-2.0');
+    writeln(psfile, '%%Creator: pldv');
+  end;
+  ps_style := 99;
+end;
+
+procedure ps_graph_off;
+begin
+end;
+
+procedure ps_draw_end;
+begin
+    writeln(psfile,'grestore');
+    writeln(psfile,'showpage');
+    close(psfile);
+end;
+
+
+procedure ps_draw_area(world_width,world_height,x,y:real;eps:boolean);
+const
+  f=2.8346;
+begin
+    if eps then begin
+      writeln(psfile,'%%BoundingBox: ',
+                                trunc(-x*f) ,' ',
+                                trunc(-y*f) ,' ',
+                                round((-x+world_width)*f),' ',
+                                round((-y+world_height)*f) );
+       writeln(psfile,'%%EndComments');
+    end;
+    writeln(psfile, 'gsave          ');
+    writeln(psfile, '2.8346 2.8346 scale ');
+    if not eps then writeln(psfile, '90 rotate           ');
+    writeln(psfile, '0.15 setlinewidth      ');
+    writeln(psfile, '1 setlinejoin       ');
+    writeln(psfile);
+    if not eps then 
+    begin
+      writeln(psfile, '/Times-Roman findfont 5 scalefont setfont ');
+      writeln(psfile, '/seite 0 def ');
+      writeln(psfile, '/str 3 string def ');
+      writeln(psfile, '/breite ',world_width:8:2,' def' );
+      writeln(psfile, '/hoehe ',world_height:8:2,' def' );
+      writeln(psfile, '/wx ',x:8:2, ' def');
+      writeln(psfile, '/wy ',y:8:2, ' def');
+      writeln(psfile,'17 wx add -9 hoehe sub wy add translate');
+      writeln(psfile,'Rahmen'); 
+    end;
+    writeln(psfile);
+end;
+
+procedure ps_newpage;
+begin
+   writeln(psfile,'gsave'); 
+   writeln(psfile,'showpage'); 
+   writeln(psfile,'grestore');
+   writeln(psfile,'Rahmen'); writeln(psfile);
+end;
+
+procedure ps_new_color(color:integer);
+begin
+  if color= default    then color := black;
+  if color= black           then writeln(psfile,'0 0 0 setrgbcolor');
+  if color= blue            then writeln(psfile,'0 0 .8 setrgbcolor');
+  if color= green           then writeln(psfile,'0 .8 0 setrgbcolor');
+  if color= cyan            then writeln(psfile,'0 .8 .8 setrgbcolor');
+  if color= red             then writeln(psfile,'.8 0 0 setrgbcolor');
+  if color= magenta         then writeln(psfile,'.8 0 .8 setrgbcolor');
+  if color= brown           then writeln(psfile,'.65 .16 .16 setrgbcolor');
+  if color= lightgray       then writeln(psfile,'.8 .8 .8 setrgbcolor');
+  if color= darkgray        then writeln(psfile,'.6 .6 .6 setrgbcolor');
+  if color= lightblue       then writeln(psfile,'0 0 1 setrgbcolor');
+  if color= lightgreen      then writeln(psfile,'0 1 0 setrgbcolor');
+  if color= lightcyan       then writeln(psfile,'0 1 1 setrgbcolor');
+  if color= lightred        then writeln(psfile,'1 0 0 setrgbcolor');
+  if color= lightmagenta    then writeln(psfile,'1 0 1 setrgbcolor');
+  if color= yellow     	    then writeln(psfile,'1 1 0 setrgbcolor');
+  if color= white           then writeln(psfile,'1 1 1 setrgbcolor');
+end;
+
+procedure ps_set_linewidth( width : real);
+begin
+  writeln(psfile,width:8:2,'  setlinewidth');
+end;
+
+procedure ps_new_linewidth(factor : real);
+const
+   default = 0.15;
+begin
+   ps_set_linewidth(default*factor);
+end;
+
+
+procedure ps_linec2d( x1,y1,x2,y2: real; style:integer);
+begin
+   ps_set_style(style);
+   writeln(psfile,x1:8:2,' ',y1:8:2,' moveto');
+   writeln(psfile,x2:8:2,' ',y2:8:2,'  lineto');
+   writeln(psfile,'stroke');
+end;
+
+procedure ps_pointc2d(x,y :real; marker:integer);
+begin
+    writeln(psfile,'/lw currentlinewidth def');
+(*    writeln(psfile,'0 setlinewidth');*)
+    writeln(psfile,'[] 0 setdash');
+    
+   (*  if marker > 3 then marker := 3; *)
+    Case marker of 
+      100 {PIXEL},	
+      3 {DOT}     : begin
+                      writeln(psfile, 'newpath');
+                      writeln(psfile, x:8:2,' ',y:8:2,' .2 0 360 arc');
+                      writeln(psfile, 'fill');
+                    end;   
+      1 {PLUS}    : begin
+                      writeln(psfile, 'newpath');
+                      writeln(psfile, x:8:2,' ',y+1:8:2,' moveto');
+                      writeln(psfile,' 0 -2 rlineto');
+                      writeln(psfile,'-1  1 rmoveto');
+                      writeln(psfile,' 2  0 rlineto');
+                      writeln(psfile,'stroke');
+                    end;  
+      10 {SMALLPLUS} : begin
+(*
+                      writeln(psfile, 'newpath');
+                      writeln(psfile, x:8:2,' ',y+1:8:2,' moveto');
+                      writeln(psfile,' 0 -1 rlineto');
+                      writeln(psfile,'-0.5  0.5 rmoveto');
+                      writeln(psfile,' 1  0 rlineto');
+                      writeln(psfile,'stroke');
+*)
+                      writeln(psfile, 'newpath');
+                      writeln(psfile, x:8:2,' ',y:8:2,' .5 0 360 arc');
+                      writeln(psfile, 'fill');
+                    end;  
+      50 {BIGPIXEL} : begin
+                      writeln(psfile, 'newpath');
+                      writeln(psfile, x:8:2,' ',y:8:2,' .35 0 360 arc');
+                      writeln(psfile, 'fill');
+                    end;  
+      2 {CROSS}   : begin
+                      writeln(psfile, 'newpath');
+                      writeln(psfile, x+0.707:8:2,' ',y+0.707:8:2,' moveto');
+                      writeln(psfile, '-1.414 -1.414 rlineto');
+                      writeln(psfile, '1.414 0 rmoveto');
+                      writeln(psfile, '-1.414  1.414 rlineto');
+                      writeln(psfile, 'stroke');
+                    end;  
+      0 {circle}  : begin
+                      writeln(psfile, 'newpath');
+                      writeln(psfile, x:8:2,' ',y:8:2,' .8 0 360 arc');
+                      writeln(psfile, 'stroke');
+                    end; 
+      end;
+    writeln(psfile,'lw setlinewidth');  
+    ps_style := 0;
+end;
+
+procedure  ps_curve2d(var p:vts2d; n1,n2,style:integer);
+var
+    i : integer;
+begin
+    ps_set_style(style);
+    writeln(psfile, 'newpath');
+    writeln(psfile, p[n1].x:8:2,' ',p[n1].y:8:2,' moveto');
+    writeln(psfile, '[');
+    for i:=n1+1 to n2 do
+        writeln(psfile,' [', p[i].x:8:2,' ',p[i].y:8:2, ']');
+    writeln(psfile, ']');
+    writeln(psfile,'{aload pop lineto} forall');
+    writeln(psfile,'stroke');    
+end;
+
+

cdg0gv-2010/include/proc_ag.pas

@@ -0,0 +1,424 @@
+{*********************************}
+{*********  P R O C _ A G  *******}
+{*********************************}
+
+function sign(a : real) : integer;
+ begin   if a<0  then  sign:= -1  else  sign:= 1;  end; {sign}
+{*************}
+function max(a,b : real) : real;
+ begin   if a>=b then  max:= a   else   max:= b;   end; {max}
+{*************}
+function min(a,b : real) : real;
+ begin   if a<=b then  min:= a   else   min:= b;   end; {min}
+{*************}
+
+procedure put2d(x,y : real; var v: vt2d);
+ begin   v.x:= x;  v.y:= y;   end;
+{*************}
+procedure put3d(x,y,z : real; var v: vt3d);
+ begin   v.x:= x;  v.y:= y;  v.z:= z;   end;
+{*************}
+
+procedure get3d(v : vt3d; var x,y,z: real);
+ begin  x:= v.x;  y:= v.y;  z:= v.z;  end;
+{*************}
+
+procedure scale2d(r : real; v: vt2d; var vs: vt2d);
+ begin   vs.x:= r*v.x;  vs.y:= r*v.y;  end;
+{*************}
+procedure scale3d(r : real; v: vt3d; var vs: vt3d);
+ begin   vs.x:= r*v.x;  vs.y:= r*v.y;  vs.z:= r*v.z;  end;
+{*************}
+
+procedure scaleco2d(r1,r2 : real; v: vt2d; var vs: vt2d);
+ begin   vs.x:= r1*v.x;  vs.y:= r2*v.y;  end;
+{*************}
+procedure scaleco3d(r1,r2,r3 : real; v: vt3d; var vs: vt3d);
+ begin   vs.x:= r1*v.x;  vs.y:= r2*v.y;  vs.z:= r3*v.z;  end;
+{*************}
+
+
+procedure sum2d(v1,v2 : vt2d; var vs : vt2d);
+ begin  vs.x:= v1.x + v2.x;  vs.y:= v1.y + v2.y; end;
+{*************}
+procedure sum3d(v1,v2 : vt3d; var vs : vt3d);
+ begin vs.x:= v1.x + v2.x;  vs.y:= v1.y + v2.y;  vs.z:= v1.z + v2.z; end;
+{*************}
+
+procedure diff2d(v1,v2 : vt2d; var vs : vt2d);
+ begin  vs.x:= v1.x - v2.x;   vs.y:= v1.y - v2.y; end;
+{*************}
+procedure diff3d(v1,v2 : vt3d; var vs : vt3d);
+ begin vs.x:= v1.x - v2.x;  vs.y:= v1.y - v2.y;  vs.z:= v1.z - v2.z; end;
+{*************}
+
+procedure lcomb2vt2d(r1: real; v1: vt2d; r2: real; v2: vt2d; var vlc : vt2d);
+ begin vlc.x:= r1*v1.x + r2*v2.x;  vlc.y:= r1*v1.y + r2*v2.y; end;
+{*************}
+procedure lcomb2vt3d(r1: real; v1: vt3d; r2: real; v2: vt3d; var vlc : vt3d);
+ begin
+ vlc.x:= r1*v1.x + r2*v2.x; vlc.y:= r1*v1.y + r2*v2.y; vlc.z:= r1*v1.z + r2*v2.z;
+ end;
+{*************}
+
+procedure lcomb3vt2d(r1: real; v1: vt2d; r2: real; v2: vt2d;
+                                         r3: real; v3: vt2d; var vlc : vt2d);
+ begin
+   vlc.x:= r1*v1.x + r2*v2.x + r3*v3.x;
+   vlc.y:= r1*v1.y + r2*v2.y + r3*v3.y;
+ end;
+{*************}
+procedure lcomb3vt3d(r1: real; v1: vt3d; r2: real; v2: vt3d;
+                                          r3: real; v3: vt3d; var vlc : vt3d);
+ begin
+   vlc.x:= r1*v1.x + r2*v2.x + r3*v3.x;
+   vlc.y:= r1*v1.y + r2*v2.y + r3*v3.y;
+   vlc.z:= r1*v1.z + r2*v2.z + r3*v3.z;
+ end;
+{*************}
+
+procedure lcomb4vt3d(r1: real; v1: vt3d; r2: real; v2: vt3d;
+                     r3: real; v3: vt3d; r4: real; v4: vt3d; var vlc : vt3d);
+ begin
+   vlc.x:= r1*v1.x + r2*v2.x + r3*v3.x + r4*v4.x;
+   vlc.y:= r1*v1.y + r2*v2.y + r3*v3.y + r4*v4.y;
+   vlc.z:= r1*v1.z + r2*v2.z + r3*v3.z + r4*v4.z;
+ end;
+{*************}
+
+function abs2d(v : vt2d) : real;
+ begin  abs2d:= abs(v.x) + abs(v.y);  end;
+{*************}
+function abs3d(v : vt3d) : real;
+ begin  abs3d:= abs(v.x) + abs(v.y) + abs(v.z);  end;
+{*************}
+
+function length2d(v : vt2d) : real;
+ begin  length2d:= sqrt( sqr(v.x) + sqr(v.y) );  end;
+{*************}
+function length3d(v : vt3d) : real;
+ begin  length3d:= sqrt( sqr(v.x) + sqr(v.y) + sqr(v.z));  end;
+{*************}
+
+procedure normalize2d(var p: vt2d);
+var c : real;
+begin c:= 1/length2d(p); p.x:= c*p.x; p.y:= c*p.y; end;
+{************}
+procedure normalize3d(var p: vt3d);
+var c : real;
+begin c:= 1/length3d(p); p.x:= c*p.x; p.y:= c*p.y; p.z:= c*p.z  end;
+{************}
+
+function scalarp2d(p1,p2 : vt2d) : real;
+ begin  scalarp2d:= p1.x*p2.x + p1.y*p2.y;  end;
+{*************}
+function scalarp3d(p1,p2 : vt3d) : real;
+ begin  scalarp3d:= p1.x*p2.x + p1.y*p2.y + p1.z*p2.z;  end;
+{*************}
+
+function distance2d(p,q : vt2d): real;
+begin  distance2d:= sqrt( sqr(p.x-q.x) + sqr(p.y-q.y) );  end;
+{*************}
+function distance3d(p,q : vt3d): real;
+begin  distance3d:= sqrt( sqr(p.x-q.x) + sqr(p.y-q.y) + sqr(p.z-q.z) );  end;
+{*************}
+function distance2d_square(p,q : vt2d) : real;
+begin distance2d_square:= sqr(p.x-q.x) + sqr(p.y-q.y);  end;
+{**************}
+function distance3d_square(p,q : vt3d) : real;
+begin distance3d_square:= sqr(p.x-q.x) + sqr(p.y-q.y) + sqr(p.z-q.z);  end;
+{**************}
+
+procedure vectorp(v1,v2 : vt3d; var vp : vt3d);
+{Berechnet das Kreuzprodukt von (x1,y1,z1) und (x2,y2,z2). }
+ begin
+   vp.x:=  v1.y*v2.z - v1.z*v2.y ;
+   vp.y:= -v1.x*v2.z + v2.x*v1.z ;
+   vp.z:=  v1.x*v2.y - v2.x*v1.y ;
+ end; {vectorp}
+{*************}
+
+function determ3d(v1,v2,v3: vt3d) : real;
+{Berechnet die Determinante einer 3x3-Matrix.}
+ begin
+   determ3d:= v1.x*v2.y*v3.z + v1.y*v2.z*v3.x + v1.z*v2.x*v3.y
+            - v1.z*v2.y*v3.x - v1.x*v2.z*v3.y - v1.y*v2.x*v3.z;
+ end; {determ3d}
+{*************}
+
+procedure rotor2d(cos_rota,sin_rota : real; p : vt2d; var pr: vt2d);
+ begin
+   pr.x:= p.x*cos_rota - p.y*sin_rota;
+   pr.y:= p.x*sin_rota + p.y*cos_rota;
+ end; {rotor2d}
+{*************}
+
+procedure rotp02d(cos_rota,sin_rota : real; p0,p : vt2d; var pr: vt2d);
+ begin
+   pr.x:= p0.x + (p.x-p0.x)*cos_rota - (p.y-p0.y)*sin_rota;
+   pr.y:= p0.y + (p.x-p0.x)*sin_rota + (p.y-p0.y)*cos_rota;
+ end; {rotor2d}
+{*************}
+
+procedure rotorz(cos_rota,sin_rota : real; p : vt3d; var pr: vt3d);
+ begin
+   pr.x:= p.x*cos_rota - p.y*sin_rota;
+   pr.y:= p.x*sin_rota + p.y*cos_rota;    pr.z:= p.z;
+ end; {rotorz}
+{*************}
+procedure rotp0z(cos_rota,sin_rota : real; p0,p: vt3d; var pr: vt3d);
+ begin
+   pr.x:= p0.x + (p.x-p0.x)*cos_rota - (p.y-p0.y)*sin_rota;
+   pr.y:= p0.y + (p.x-p0.x)*sin_rota + (p.y-p0.y)*cos_rota;    pr.z:= p.z;
+ end; {rotp0z}
+{*************}
+procedure rotp0x(cos_rota,sin_rota : real; p0,p: vt3d; var pr: vt3d);
+{Rotation um eine zur x-Achse parallele Achse durch p0}
+ begin
+   pr.y:= p0.y + (p.y-p0.y)*cos_rota - (p.z-p0.z)*sin_rota;
+   pr.z:= p0.z + (p.y-p0.y)*sin_rota + (p.z-p0.z)*cos_rota;    pr.x:= p.x;
+ end; {rotp0x}
+{*************}
+procedure rotp0y(cos_rota,sin_rota : real; p0,p: vt3d; var pr: vt3d);
+{Rotation um eine zur y-Achse parallele Achse durch p0}
+ begin
+   pr.x:= p0.x + (p.x-p0.x)*cos_rota - (p.z-p0.z)*sin_rota;
+   pr.z:= p0.z + (p.x-p0.x)*sin_rota + (p.z-p0.z)*cos_rota;    pr.y:= p.y;
+ end; {rotp0y}
+{*************}
+
+procedure change1d(var a,b: real);
+var aa: real;
+begin  aa:= a;  a:= b;  b:= aa;  end;
+{*********}
+procedure change2d(var v1,v2: vt2d);
+var vv: vt2d;
+begin  vv:= v1;  v1:= v2;  v2:= vv;  end;
+{**********}
+procedure change3d(var v1,v2: vt3d);
+var vv: vt3d;
+begin  vv:= v1;  v1:= v2;  v2:= vv;  end;
+{**********}
+
+{polar_angles fuer altes p2c: (SUSE 6.4)}
+function polar_angle(x,y : real) : real;
+ var w : real;
+ begin
+    if (x=0) and (y=0) then   w:= 0
+     else
+      begin
+        if abs(y)<=abs(x) then
+           begin
+             w:= arctan(y/x);
+             if x<0 then  w:=pi+w
+                    else  if (y<0) and( w<>0) then  w:= pi2+w;
+           end
+          else
+           begin
+             w:= pih-arctan(x/y);
+             if y<0 then w:= pi+w;
+           end; {if}
+      end;  {if}
+    polar_angle:= w;
+ end;  { polar_angle }
+{******************}
+
+procedure equation_degree1(a,b :real; var x : real ; var ns : integer);
+ {Gleichung 1.Grades: a*x + b = 0 }
+ begin
+  if abs(a)=0 {<eps8} then   ns:= -1
+                 else   begin x:=-b/a;   ns:= 1;  end;
+ end; { equation_degree1 }
+{***************}
+
+procedure equation_degree2(a,b,c:real; var x1,x2:real ;var ns:integer);
+{Berechnet die REELLEN Loesungen einer Gleichnung 2.Grades: a*x*x+b*x+c = 0.}
+ var dis,wu2,xx1 : real;
+ begin
+   ns:=0;
+   if abs(a)=0 {<eps8} then
+        equation_degree1(b,c,x1,ns)
+      else
+        begin
+        dis:= b*b-4*a*c;
+        if (dis<0) and (dis>-eps6*abs(b)) then  dis:=0;
+        if dis<0 then   ns:= 0
+                 else
+                   begin
+                   if abs(dis)<eps8 then
+                      begin  ns:= 1;   x1:= -b/(2*a);  x2:= x1; end
+                     else
+                      begin
+                      ns:= 2;                wu2:= sqrt(dis);
+                      x1:= (-b+wu2)/(2*a);    x2:= (-b-wu2)/(2*a);
+                      end;
+                   end; {dis>=0}
+        end; {a<>0}
+{Umordnen nach Groesse:}
+   if ns=2 then  if x2<x1 then  change1d(x1,x2);
+ end;  { equation_degree2 }
+{*****************}
+
+procedure is_line_line(a1,b1,c1, a2,b2,c2 : real; var xs,ys : real;
+                                                  var nis : integer);
+ {Schnittpunkt (xs,ys) (ns=1) der Geraden a1*x+b1*y=c1, a2*x+b2*y=c2 .}
+ var det : real;
+ begin
+   det:= a1*b2 - a2*b1;
+   if abs(det)<eps8 then    nis:= 0
+       else
+         begin
+         nis:= 1;  xs:= (c1*b2 - c2*b1)/det;   ys:= (a1*c2 - a2*c1)/det;
+         end;
+ end;   { is_line_line }
+{************}
+
+procedure is_unitcircle_line(a,b,c: real; var x1,y1,x2,y2 : real;
+                                                  var nis : integer);
+ {Schnitt Kreis-Gerade: sqr(x)+sqr(y)=1, a*x+b*y=c,
+  Schnittpkte: (x1,y1),(x2,y2).Es ist x1<=x2, nis Anzahl der Schnittpunkte}
+ var ab2,wu,dis : real;
+ begin
+   nis:= 0;   ab2:= a*a + b*b;   dis:= ab2 - c*c;
+   if dis>=0 then
+      begin
+      nis:= 2;
+      wu:= sqrt(dis);           if abs(wu)<eps8 then nis:= 1;
+      x1:= (a*c-b*wu)/ab2;      y1:= (b*c+a*wu)/ab2;
+      x2:= (a*c+b*wu)/ab2;      y2:= (b*c-a*wu)/ab2;
+      if x2<x1 then begin  change1d(x1,x2);  change1d(y1,y2);  end;
+      end;
+ end;   { is_unitcircle_line }
+{************}
+
+procedure is_circle_line(xm,ym,r, a,b,c: real; var x1,y1,x2,y2 : real;
+                                               var nis : integer);
+ {Schnitt Kreis-Gerade: sqr(x-xm)+sqr(y-ym)=r*r, a*x+b*y=c,
+  Schnittpkte: (x1,y1),(x2,y2).Es ist x1<=x2, nis Anzahl der Schnittpunkte}
+ var ab2,wu,dis,cc : real;
+ begin
+   nis := 0;
+   ab2:= a*a + b*b;   cc:=c - a*xm - b*ym;   dis:= r*r*ab2 - cc*cc;
+   if dis>=0 then
+      begin
+      nis:= 2;
+      wu:= sqrt(dis);           if abs(wu)<eps8 then nis:= 1;
+      x1:= xm + (a*cc-b*wu)/ab2;   y1:= ym + (b*cc+a*wu)/ab2;
+      x2:= xm + (a*cc+b*wu)/ab2;   y2:= ym + (b*cc-a*wu)/ab2;
+      if x2<x1 then begin  change1d(x1,x2);  change1d(y1,y2);  end;
+      end;
+ end;   { is_circle_line }
+{************}
+
+procedure is_circle_circle(xm1,ym1,r1,xm2,ym2,r2: real;
+                                 var x1,y1,x2,y2: real; var nis: integer);
+{Schnitt Kreis-Kreis. Es ist x1<x2. nis = Anzahl der Schnittpunkte }
+ var a,b,c : real;
+ begin
+   a:= 2*(xm2 - xm1);        b:= 2*(ym2 - ym1);
+   c:= r1*r1 - sqr(xm1) -sqr(ym1) - r2*r2 + sqr(xm2) + sqr(ym2);
+   if (abs(a)+abs(b)<eps8) then   nis:= 0
+      else   is_circle_line(xm1,ym1,r1,a,b,c,x1,y1,x2,y2,nis);
+ end;   { is_circle_circle }
+{************}
+
+function pt_before_plane(p,nv: vt3d; d: real) : boolean;
+{...stellt fest, ob der Punkt p "vor" der Ebene  nv*x-d=0 liegt.}
+ begin  if (scalarp3d(p,nv)-d) >=0 then pt_before_plane:= true
+                                   else pt_before_plane:= false;  end;
+{*************}
+
+procedure plane_equ(p1,p2,p3 : vt3d; var nv : vt3d; var d: real;
+                                                      var error : boolean);
+{Berechnet die Gleichung nv*x=d der Ebene durch die Punkte p1,p2,p3.
+ error=true: die Punkte spannen keine Ebene auf.}
+ var p21,p31 : vt3d;
+ begin
+   diff3d(p2,p1,p21);      diff3d(p3,p1,p31);
+   vectorp(p21,p31,nv);    d:= scalarp3d(nv,p1);
+   if abs3d(nv)<eps8 then error:= true else error:= false;
+ end;   { plane_coeff }
+{*************}
+
+procedure is_line_plane(p,rv,nv: vt3d; d : real; var pis : vt3d;
+                                                 var nis : integer);
+{Schnitt Gerade-Ebene. Gerade: Punkt p, Richtung r. Ebene:  nv*x = d .
+ nis=0: kein Schnitt ,nis=1: Schnittpunkt, nis=2: Gerade liegt in der Ebene.}
+ var t,sp,pd : real;
+ begin
+   sp:= scalarp3d(nv,rv);   pd:= scalarp3d(nv,p) - d;
+   if abs(sp)<eps8 then
+        begin  nis:= 0;  if abs(pd)<eps8 then  nis:= 2;  end
+       else
+        begin
+        t:= -pd/sp;   lcomb2vt3d(1,p ,t,rv ,pis);    nis:= 1;
+        end;
+ end;   { is_line_plane }
+{*************}
+
+procedure is_3_planes(nv1: vt3d; d1: real; nv2: vt3d; d2: real;
+                      nv3: vt3d; d3: real;
+                      var pis : vt3d; var error: boolean);
+{Schnitt der Ebenen nv1*x=d1, nv2*x=d2, nv3*x=d3.
+ error= true: Schnitt besteht nicht aus einem Punkt.}
+ var det,dd1,dd2,dd3 : real;  n12,n23,n31 : vt3d;
+ begin
+   vectorp(nv1,nv2, n12);  vectorp(nv2,nv3, n23);  vectorp(nv3,nv1, n31);
+   det:= scalarp3d(nv1,n23);
+   if abs(det)<eps8 then  error:= true
+      else
+       begin
+       dd1:= d1/det;  dd2:= d2/det;  dd3:= d3/det;
+       lcomb3vt3d(dd1,n23, dd2,n31, dd3,n12 ,pis);
+       end;
+ end;   { is_3_planes }
+{*************}
+
+procedure is_plane_plane(nv1: vt3d; d1: real; nv2: vt3d; d2: real;
+                             var p,rv : vt3d; var error: boolean);
+{Schnitt der Ebenen nv1*x=d1, nv2*x=d2. Schnittgerade: x = p +  t*rv .
+ error= true: Schnitt besteht nicht aus einem Punkt.}
+ var det,c11,c22,c12,s1,s2: real;
+ begin
+   c11:= scalarp3d(nv1,nv1);     c22:= scalarp3d(nv2,nv2);
+   c12:= scalarp3d(nv1,nv2);     det:= c11*c22-sqr(c12);
+   if abs(det)=0 {<eps8} then  error:= true
+      else
+       begin
+       s1:= (d1*c22-d2*c12)/det;        s2:= (d2*c11-d1*c12)/det;
+       lcomb2vt3d(s1,nv1, s2,nv2, p);   vectorp(nv1,nv2, rv);
+       error:= false;
+       end;
+ end;   { is_plane_plane }
+{***********}
+
+procedure ptco_plane3d(p0,v1,v2,p: vt3d; var xi,eta: real; var error: boolean);
+ {v1,v2 sind linear unabhaengig, p-p0 linear abhaengig von v1,v2.
+ Es werden Zahlen  xi,eta berechnet mit  p = p0 + xi*v1 + eta*v2. }
+ var det,s11,s12,s22,s13,s23 : real;   v3 : vt3d;
+ begin
+   diff3d(p,p0, v3);
+   s11:= scalarp3d(v1,v1);   s12:= scalarp3d(v1,v2);   s13:= scalarp3d(v1,v3);
+                             s22:= scalarp3d(v2,v2);   s23:= scalarp3d(v2,v3);
+   det := s11*s22 - s12*s12;
+   if abs(det)=0 {<eps8} then  error:= true
+      else
+        begin
+        error:= false;
+        xi := (s13*s22 - s23*s12)/det;
+        eta:= (s11*s23 - s13*s12)/det;
+        end;
+ end;   { ptco_plane3d }
+{*****************}
+
+
+procedure newcoordinates3d(p,b0,b1,b2,b3: vt3d; var pnew: vt3d);
+{Berechnet die Koordinaten von p bzgl. der Basis b1,b2,b3 mit Nullpkt. b0.}
+var det : real;   p0: vt3d;
+begin
+  diff3d(p,b0, p0);              det:= determ3d(b1,b2,b3);
+  pnew.x:= determ3d(p0,b2,b3)/det;
+  pnew.y:= determ3d(b1,p0,b3)/det;
+  pnew.z:= determ3d(b1,b2,p0)/det;
+end;  { newcoordinates3d }  
+{***}

cdg0gv-2010/include/proc_pp.pas

@@ -0,0 +1,196 @@
+{**********************}
+{***  P R O C  P P  ***}
+{**********************}
+
+procedure init_parallel_projection;
+ begin
+   writeln('*** PARALLEL-PROJEKTION ***');
+   writeln;
+   writeln('Projektionswinkel u, v ? (in Grad)');  readln(u_angle,v_angle);
+   rad_u:= u_angle*pi/180  ;  rad_v:= v_angle*pi/180  ;
+   sin_u:= sin(rad_u)   ;           cos_u:= cos(rad_u)   ;
+   sin_v:= sin(rad_v)   ;           cos_v:= cos(rad_v)   ;
+{Normalen-Vektor der Bildebene:}
+   n0vt.x:= cos_u*cos_v;    n0vt.y:= sin_u*cos_v;    n0vt.z:= sin_v;
+ end; { init_parallel_projection }
+{**************}
+
+procedure pp_vt3d_vt2d(p : vt3d; var pp: vt2d);
+ {Bildvektor eines Punktes}
+ begin
+   pp.x:= -p.x*sin_u + p.y*cos_u;
+   pp.y:=-(p.x*cos_u + p.y*sin_u)*sin_v + p.z*cos_v;
+ end; {pp_vt3d_vt2d}
+{*************}
+
+procedure pp_point(p: vt3d; style: integer);
+ {... markiert einen projizierten Punkt}
+ var pp : vt2d;
+ begin
+   pp_vt3d_vt2d(p,pp);   point2d(pp,style);
+ end; {pp_point}
+{*************}
+
+procedure pp_line(p1,p2 : vt3d ; style : integer);
+ {projiziert die Strecke p1,p2 }
+ var pp1,pp2 : vt2d;
+ begin
+   pp_vt3d_vt2d(p1,pp1);       pp_vt3d_vt2d(p2,pp2);
+   line2d(pp1,pp2,style);
+ end; {pp_line}
+{*************}
+
+procedure pp_arrow(p1,p2 : vt3d; style : integer);
+ {projiziert einen Pfeil}
+ var pp1,pp2: vt2d;
+ begin
+   pp_vt3d_vt2d(p1,pp1);   pp_vt3d_vt2d(p2,pp2);
+   arrow2d(pp1,pp2,style);
+ end; {pp_axes}
+{*************}
+
+procedure pp_axes(al : real);
+ {projiziert die Koordinatenachsen, al: Achsenlaenge}
+ var p0,p1,p2,p3 : vt3d;
+ begin
+   put3d(0,0,0,p0);  put3d(al,0,0,p1);
+   put3d(0,al,0,p2); put3d(0,0,al,p3);
+   pp_arrow(p0,p1, 2);   pp_arrow(p0,p2, 2);
+   pp_arrow(p0,p3, 2);
+ end; {pp_axes}
+{*************}
+
+procedure pp_vts3d_vts2d(var p: vts3d; n1,n2 : integer; var pp : vts2d);
+ {Berechnet die Bildvektoren pp einer Punktreihe p.}
+ var i : integer;
+  begin
+    for i:= n1 to n2 do
+      begin
+        pp[i].x:= - p[i].x*sin_u + p[i].y*cos_u;
+        pp[i].y:= -(p[i].x*cos_u + p[i].y*sin_u)*sin_v + p[i].z*cos_v;
+      end;
+  end; {pp_vts3d_vts2d}
+{*************}
+
+procedure pp_curve(var p: vts3d; n1,n2,style : integer);
+ {Projiziert das 3d-Polygon p[n1]...p[n2]}
+ var pp : vts2d;
+ begin
+   pp_vts3d_vts2d(p,n1,n2,pp);
+   curve2d(pp,n1,n2,style);
+ end; {pp_curve}
+{*************}
+
+procedure pp_curve_vis(var p : vts3d; n1,n2,style : integer; visible: b_array);
+ {projiziert ein 3d-Polygon. Es werden je zwei benachnarte "visible"
+ Punkte verbunden. style=10: Rest wird gestrichelt.}
+ var pp : vts2d;
+ begin
+   pp_vts3d_vts2d(p,n1,n2,pp);
+   curve2d_vis(pp,n1,n2,style,visible);
+ end; {pp_curve_vis}
+{*************}
+
+procedure pp_line_before_plane(p1,p2,nv: vt3d; d : real; side,style: integer);
+ {Zeichnet ,falls style=0, den Teil der Strecke p1, p2, f<EFBFBD>r den
+ side*(nv*x - d) >= 0 ist und strichelt den Rest, falls style = 10 ist. }
+ var dis1,dis2,t : real;    p3 : vt3d;
+ begin
+   dis1:= scalarp3d(nv,p1) - d;   dis2:= scalarp3d(nv,p2) - d;
+   if side<0 then begin  dis1:= -dis1; dis2:= -dis2;  end;
+   if (dis1>=0) and (dis2>=0) then  pp_line(p1,p2,0)
+     else
+        if (dis1<0) and (dis2<0) then
+            begin
+            if style=10 then      pp_line(p1,p2,1);
+            end
+          else
+            begin
+            t := -dis1/(dis2-dis1);
+            lcomb2vt3d(1-t,p1, t,p2, p3);    {Schnittpunkt}
+            if dis1>=0    then           {p1 vor der Ebene}
+               begin
+               pp_line(p1,p3,0);
+               if style=10 then pp_line(p3,p2,1);
+               end
+              else
+               begin
+               pp_line(p3,p2,0);
+               if style=10 then pp_line(p1,p3,1);
+               end;
+            end;
+  end; { pp_line_before_plane }
+{**************}
+
+procedure pp_curve_before_plane(var p: vts3d ; n1,n2: integer; nv: vt3d;
+                                         d: real; side,style : integer);
+ {Projiziert eine Kurve (Polygonzug) "vor" der Ebene nv*x = d und strichelt
+ den Rest, falls style=10 .}
+ var i,k,nis : integer;  nvv,p_i,p_i1,rvt,pis : vt3d;    visi,visi1 : boolean;
+     pp : vts3d;   vis : b_array;       dd : real;
+ begin
+   if side<0 then begin scale3d(-1,nv,nvv);  dd:= -d;  end
+             else begin nvv:= nv;  dd:= d; end;
+   k:= n1-1;
+   for i:= n1 to n2 do
+       begin
+       p_i:= p[i]; k:= k+1;
+       if scalarp3d(p_i,nvv)-dd >= 0 then visi:= true
+                                    else visi:= false;
+       pp[k]:= p_i;  vis[k]:= visi;
+       if (i>n1) and (visi<>visi1) then
+             begin
+             diff3d(p_i,p_i1, rvt);
+             is_line_plane(p_i1,rvt, nvv,dd, pis,nis);
+             k:= k+1;  pp[k]:= pis;  vis[k]:= true;
+             end;
+       p_i1:= p_i;  visi1:= visi;
+       end;    { for i }
+  pp_curve_vis(pp,n1,k,style,vis);
+ end; { pp_curve_before_plane }
+{*************}
+(*
+
+procedure pp_curve_before_plane(var p: vts3d ; n1,n2: integer; nv: vt3d; d: real;
+                                                       side,style : integer);
+{Projiziert eine Kurve (Polygonzug) "vor" der Ebene nv*x = d und strichelt
+ den Rest, falls style=10 .}
+ var i : integer;
+ begin
+   for i:= n1 to n2-1  do
+       pp_line_before_plane(p[i],p[i+1],nv,d,side,style);
+ end; { pp_curve_before_plane }
+{*************}
+
+procedure pp_curve_before_plane(var p: vts3d ; n1,n2: integer; nv: vt3d;
+                                         d: real; side,style : integer);
+ {Projiziert eine Kurve (Polygonzug) "vor" der Ebene nv*x = d und strichelt
+ den Rest, falls style=10 .}
+ var i : integer;  nvv : vt3d;  pp1,pp2 : vt2d;
+     vis1,vis2 : boolean;            dd : real;
+ begin
+   if side<0 then begin scale3d(-1,nv,nvv);  dd:= -d;  end
+             else begin nvv:= nv;  dd:= d; end;
+   for i:= n1 to n2 do
+       begin
+       if scalarp3d(p[i],nvv)-dd >= 0 then
+          begin
+          vis2:=true;    pp_vt3d_vt2d(p[i],pp2);
+          end
+         else
+          begin
+          vis2:= false;  if style=10 then pp_vt3d_vt2d(p[i],pp2);
+          end;
+       if i>n1 then
+          begin
+          if vis1 and vis2 then line2d(pp1,pp2,style);
+          if (style=10) and (not vis1) and (not vis2) then line2d(pp1,pp2,1);
+          if (vis1 and not vis2) or (not vis1 and vis2)
+             then pp_line_before_plane(p[i-1],p[i],nv,d,side,style);
+          end;
+       pp1:= pp2;  vis1:= vis2;
+       end;
+ end; { pp_curve_before_plane }
+{*************}
+
+*)

cdg0gv-2010/include/proc_zp.pas

@@ -0,0 +1,171 @@
+{*******************************}
+{*****  P R O C _ C P      *****}
+{*****  Zentralprojektion  *****}
+{*******************************}
+
+procedure init_central_projection;
+ begin
+   writeln('*** CENTRAL-PROJECTION ***');
+   writeln;
+   writeln('mainpoint ?');  readln(mainpt.x,mainpt.y,mainpt.z);
+   writeln('distance ?');     readln(distance);
+   writeln('angles u, v ? (in degree)');  readln(u_angle,v_angle);
+   rad_u:= u_angle*pi/180;        rad_v:= v_angle*pi/180;
+   sin_u:= sin(rad_u);            cos_u:= cos(rad_u);
+   sin_v:= sin(rad_v);            cos_v:= cos(rad_v);
+{base e1,e2 and normal n0 of image plane:}
+   e1vt.x:= -sin_u;         e1vt.y:= cos_u;          e1vt.z:= 0;
+   e2vt.x:= -cos_u*sin_v;   e2vt.y:=-sin_u*sin_v;    e2vt.z:= cos_v;
+   n0vt.x:= cos_u*cos_v;    n0vt.y:= sin_u*cos_v;    n0vt.z:= sin_v;
+{centre:}
+   lcomb2vt3d(1,mainpt, distance,n0vt,  centre);
+ end; { init_central_projection }
+{**************}
+
+procedure init_centralparallel_projection(ind : integer);
+ begin
+   if ind=1 then
+             begin
+             writeln('*** CENTRAL-projection ***');    
+             writeln('mainpoint ?');  readln(mainpt.x,mainpt.y,mainpt.z);
+             writeln('distance ?');     readln(distance);
+             end 
+            else
+             begin
+             writeln('*** PARALLEL-projection ***');
+             mainpt:= null3d;   distance:= 1000000000;
+             end;
+   writeln('angles u, v ? (in degree)');  readln(u_angle,v_angle);
+   rad_u:= u_angle*pi/180;        rad_v:= v_angle*pi/180;
+   sin_u:= sin(rad_u);            cos_u:= cos(rad_u);
+   sin_v:= sin(rad_v);            cos_v:= cos(rad_v);
+{base e1,e2 and normal n0 of image plane:}
+   e1vt.x:= -sin_u;         e1vt.y:= cos_u;          e1vt.z:= 0;
+   e2vt.x:= -cos_u*sin_v;   e2vt.y:=-sin_u*sin_v;    e2vt.z:= cos_v;
+   n0vt.x:= cos_u*cos_v;    n0vt.y:= sin_u*cos_v;    n0vt.z:= sin_v;
+{centre:}
+   lcomb2vt3d(1,mainpt, distance,n0vt,  centre);
+ end; { init_centralparallel_projection }
+{**************}
+
+procedure transf_to_e1e2n0_base(p : vt3d; var pm : vt3d);
+{coordinates in system "mainpoint, e1,e2,n0".}
+ var pd : vt3d;
+ begin
+   diff3d(p,mainpt, pd);
+   pm.x:= scalarp3d(pd,e1vt);
+   pm.y:= scalarp3d(pd,e2vt);
+   pm.z:= scalarp3d(pd,n0vt);
+ end;  { transf_to_e1e2n0_base }
+{***************}
+
+procedure cp_vt3d_vt2d(p: vt3d; var pp : vt2d);
+{central projetion (coordinates) of a point}
+ var xe,ye,ze,cc : real;   pm : vt3d;
+ begin
+   diff3d(p,mainpt, pm);
+   xe:= scalarp3d(pm,e1vt);   { coordinates von p in system: }
+   ye:= scalarp3d(pm,e2vt);   { origin = mainpoint           }
+   ze:= scalarp3d(pm,n0vt);   { and base vectors  e1,e2,n0   }
+   cc:= 1-ze/distance;
+   if cc>eps6 then begin  pp.x:= xe/cc; pp.y:= ye/cc; end      { projection }
+              else
+	      writeln('central proj.: point not in front of vanishing plane !!');
+ end;  { cp_vt3d_vt2d }
+{**************}
+
+procedure cp_point(p: vt3d; style: integer);
+ var pp : vt2d;
+ begin
+   cp_vt3d_vt2d(p,pp);   point2d(pp,style);
+ end;  { cp_point }
+{**************}
+
+procedure cp_line(p1,p2 : vt3d ; style : integer);
+ var pp1,pp2 : vt2d;
+ begin
+   cp_vt3d_vt2d(p1,pp1);   cp_vt3d_vt2d(p2,pp2);   line2d(pp1,pp2,style);
+ end; {cp_line}
+{**************}
+
+procedure cp_arrow(p1,p2 : vt3d; style : integer);
+ var pp1,pp2: vt2d;
+ begin
+   cp_vt3d_vt2d(p1,pp1);   cp_vt3d_vt2d(p2,pp2);   arrow2d(pp1,pp2,style);
+ end;  { cp_arrow }
+{**************}
+
+procedure cp_axes(al : real);
+ var p0,p1,p2,p3 : vt3d;
+ begin
+   put3d(0,0,0,p0);  put3d(al,0,0,p1);
+   put3d(0,al,0,p2); put3d(0,0,al,p3);
+   cp_arrow(p0,p1, 2);   cp_arrow(p0,p2, 2);   cp_arrow(p0,p3, 2);
+ end;  { cp_axes }
+{**************}
+
+procedure cp_vts3d_vts2d(var p: vts3d; n1,n2 : integer; var pp : vts2d);
+ var i : integer;
+ begin
+   for i:= n1 to n2 do  cp_vt3d_vt2d(p[i],pp[i]);
+ end;  { cp_vts3d_vts2d }
+{*************}
+
+procedure cp_curve(var p: vts3d; n1,n2,style : integer);
+ var pp : vts2d;
+ begin
+   cp_vts3d_vts2d(p,n1,n2,pp);   curve2d(pp,n1,n2,style);
+ end;  { cp_curve }
+{*************}
+
+procedure cp_curve_vis(var p : vts3d; n1,n2,style : integer; visible: b_array );
+{connects neighbored visible points.}
+ var pp : vts2d;
+ begin
+   cp_vts3d_vts2d(p,n1,n2,pp);   curve2d_vis(pp,n1,n2,style,visible);
+ end;  { cp_curve_vis }
+{*************}
+
+procedure cp_line_before_plane(p1,p2,nv: vt3d; d : real; side,style: integer);
+{style=0: draws visible part of the line p1p2 (visible:  side*(nv*x - d) >= 0).
+ Dash-dots the invisible part if style = 10 .}
+var dis1,dis2,t : real;    p3 : vt3d;
+ begin
+   dis1:= scalarp3d(nv,p1) - d;   dis2:= scalarp3d(nv,p2) - d;
+   if side<0 then begin  dis1:= -dis1; dis2:= -dis2;  end;
+   if (dis1>=0) and (dis2>=0) then  cp_line(p1,p2,0)
+     else
+        if (dis1<0) and (dis2<0) then
+           begin
+           if style=10 then      cp_line(p1,p2,1);
+           end
+          else
+           begin
+           t := -dis1/(dis2-dis1);
+           lcomb2vt3d(1-t,p1, t,p2, p3);    {Schnittpunkt}
+           if dis1>=0    then           {p1 vor der Ebene}
+              begin
+              cp_line(p1,p3,0);
+              if style=10 then cp_line(p3,p2,1);
+              end
+             else
+              begin
+              cp_line(p3,p2,0);
+              if style=10 then cp_line(p1,p3,1);
+              end;
+            end;
+  end;  { cp_line_before_plane }
+{**************}
+
+procedure cp_curve_before_plane(var p: vts3d ; n1,n2: integer; nv: vt3d; d: real;
+                                                       side,style : integer);
+{see cp_line_before_plane}
+ var i : integer;
+ begin
+   for i:= n1 to n2-1  do
+       cp_line_before_plane(p[i],p[i+1],nv,d,side,style);
+ end;  { cp_curve_before_plane }
+{*************}
+
+
+

cdg0gv-2010/include/proc_zpo.pas

@@ -0,0 +1,668 @@
+{*****************}
+{***  PROC_ZPO ***}
+{*****************}
+
+procedure aux_polyhedron;
+{ne, nf : Anzahl der Punkte,Kanten,Flaechen
+ face[i].npf : Anzahl der Punkte (Kanten) der i-ten Flaeche
+ edge.ep1, ep2 : Anfangs- bzw. Endpunkt der k-ten Kante
+ face[i].fp[k] : k-ter Punkt der i-ten Flaeche (positiv orientiert!!!)
+ face[i].fe[k] : k-te Kante der i-ten Flaeche
+ ! Dieses UP berechnet aus nf,face[i].fp und face[i].npf:
+ ne, face[i].fe[k], edge[i].ep1 und edge[i].ep2 !.}
+ type sedge  = record  p1,p2,nr: integer;  end;
+      sedges = array[1..nemax] of sedge;
+ var i,k,i1,i2,i_edge : integer;
+     new_edge : sedge;
+     sorted_edges: sedges;
+     new_is_old,error: boolean;
+{***}
+procedure search_insert_edge(var sorted_edges: sedges; new_edge: sedge;
+                             var i_edge: integer; var new_is_old: boolean);
+{Prueft, ob new_edge in sorted_edges enthalten ist und fuegt new_edge
+ gegebenenfalls ein.}
+var found: boolean;  left,right,middle: integer;   i : integer;
+{**}
+function edge1_smaller_edge2(var edge1,edge2 : sedge) : boolean;
+begin
+  if (edge1.p1<edge2.p1) or ((edge1.p1=edge2.p1) and (edge1.p2<edge2.p2))
+     then  edge1_smaller_edge2:= true  else  edge1_smaller_edge2:= false;
+end;  { edge1_smaller_edge2 }
+{**}
+begin  { search_insert_edge }
+  left:= 1;  right:= ne;
+  if ne=0 then begin found:= false; middle:= 0; end
+     else           { binaeres Suchen }
+      repeat
+      middle:= (left+right) div 2;
+      if edge1_smaller_edge2(new_edge,sorted_edges[middle])
+         then  right:= middle-1  else  left:= middle+1;
+      with sorted_edges[middle] do if (new_edge.p1=p1) and (new_edge.p2=p2)
+                                      then found:=true else found:= false;
+      until found or (left>right);
+  if found then begin i_edge:= middle; new_is_old:= true; end
+     else
+      begin
+      new_is_old:= false; if left>middle then middle:= middle+1;
+      for i:= ne downto middle do sorted_edges[i+1]:= sorted_edges[i];
+      sorted_edges[middle]:= new_edge;
+      sorted_edges[middle].nr:= ne+1;
+      end;  { if }
+end;   { search_insert_edge }
+{***}
+ begin   { aux_polyhedron }
+   ne:= 0;
+   for i:= 1 to nf do
+       with face[i] do
+       begin
+       i1:= fp[1];
+       for k:= 1 to npf do
+           begin
+           if k<npf then i2:= fp[k+1] else i2:= fp[1];
+           new_edge.p1:= i1;  new_edge.p2:= i2;
+           with new_edge do if p1>p2 then begin p1:= i2; p2:= i1; end;
+           search_insert_edge(sorted_edges,new_edge, i_edge,new_is_old);
+           if new_is_old then fe[k]:= sorted_edges[i_edge].nr
+              else
+               begin
+               ne:= ne+1; fe[k]:= ne;
+               with edge[ne] do begin  ep1:= i1; ep2:= i2; vis:= true end;
+               end;  { if }
+           i1:= i2;
+           end;  { for k }
+       nef:= npf;
+       plane_equ(p[fp[1]],p[fp[2]],p[fp[3]], nv,d,error);
+       end; { with face }
+ end; { aux_lines_before_faces }
+{*************}
+
+procedure cp_vts3d_vts2d_spez(var p: vts3d; n1,n2: integer;
+                              var pp: vts2d; var pdist: r_array);
+{Zentralprojektion (Koordinaten) einer Punktreihe.
+ pdist[i]: Distanz des Punktes p[i] von der Bildtafel. }
+ var xe,ye,ze,cc : real;   pm : vt3d;   i : integer;
+ begin
+   for i:= n1 to n2 do
+       begin
+   diff3d(p[i],mainpt, pm);
+   xe:= scalarp3d(pm,e1vt);   { Koordinaten von p bzgl. dem Koord.-System: }
+   ye:= scalarp3d(pm,e2vt);   { Nullpunkt = Hauptpunkt                     }
+   ze:= scalarp3d(pm,n0vt);   { und Basis e1,e2,n0                         }
+   cc:= 1-ze/distance;   pdist[i]:= ze;  { Abstand zur Bildebene }
+   if cc>eps6 then begin  pp[i].x:= xe/cc; pp[i].y:= ye/cc; end
+              else
+              writeln('Punkt liegt in oder hinter der Verschwindungsebene !!');
+       end; { for }
+ end; { cp_vts3d_vts2d_spez }
+{**************}
+
+procedure aux_quadrangle(n1,n2,np0,ne0,nf0: integer);
+ var  ne1,i,k : integer;  error : boolean;
+{***}
+function f(i,k :integer) : integer;  begin f := (k-1)*(n1-1)+i+nf0;  end;
+function pt(i,k:integer) : integer;  begin pt:= (k-1)*n1+i+np0; end;
+function e1(i,k:integer) : integer;  begin e1:= (k-1)*(n1-1)+i+ne0; end;
+function e2(i,k:integer) : integer;  begin e2:= ne1+(i-1)*(n2-1)+k+ne0; end;
+{****}
+ begin
+   ne1:= (n1-1)*n2;
+   for k:= 1 to n2 do
+       begin
+       for i:= 1 to n1 do
+           begin
+           if (i<n1) and (k<n2) then
+           with face[f(i,k)] do
+                begin
+                fp[1]:= pt(i,k);      fp[2]:= pt(i+1,k);
+                fp[3]:= pt(i+1,k+1);  fp[4]:= pt(i,k+1);
+                fe[1]:= e1(i,k);      fe[2]:= e2(i+1,k);
+                fe[3]:= e1(i,k+1);    fe[4]:= e2(i,k);
+                npf:= 4; nef:= 4;
+                end;  { with }
+           if i<n1 then
+           with edge[e1(i,k)] do
+                begin  ep1:= pt(i,k);  ep2:= pt(i+1,k);  vis:= true; end;
+           if k<n2 then
+           with edge[e2(i,k)] do
+                begin  ep1:= pt(i,k);  ep2:= pt(i,k+1);  vis:= true; end;
+           end;  { for i }
+       end;  { for k }
+   np:= np0+n1*n2;  ne:= ne0+n1*(n2-1)+(n1-1)*n2;  nf:= nf0+(n1-1)*(n2-1);
+   for i:= nf0+1 to nf do
+       with face[i] do plane_equ(p[fp[1]],p[fp[2]],p[fp[3]], nv,d,error);
+ end;  { aux_quadrangle }
+{************}
+
+procedure aux_quadrangle_triang(n1,n2: integer;
+                                show_triangles: boolean);
+{Berechnet ne,nf und fuer jede Flaeche fp,fe und jede Kante ep1,ep2.}
+ var  ne1,i,k : integer;
+{***}
+function f(i,k,u:integer): integer;  begin f := (k-1)*(n1-1)+i + u*nf;  end;
+function pt(i,k:integer) : integer;  begin pt:= (k-1)*n1+i; end;
+function e1(i,k:integer) : integer;  begin e1:= (k-1)*(n1-1)+i; end;
+function e2(i,k:integer) : integer;  begin e2:= ne1+(i-1)*(n2-1)+k; end;
+function e3(i,k:integer) : integer;  begin e3:= ne+(k-1)*(n1-1)+i; end;
+{****}
+ begin
+   ne:= n1*(n2-1)+(n1-1)*n2;  nf:= (n1-1)*(n2-1);
+   ne1:= (n1-1)*n2;
+   for k:= 1 to n2-1 do
+       begin
+      for i:= 1 to n1-1 do
+          begin
+          with face[f(i,k,0)] do
+               begin
+               fp[1]:= pt(i,k);   fp[2]:= pt(i+1,k);   fp[3]:= pt(i+1,k+1);
+               fe[1]:= e1(i,k);   fe[2]:= e2(i+1,k);   fe[3]:= e3(i,k);
+               end;  { with }
+          with face[f(i,k,1)] do
+               begin
+               fp[1]:= pt(i,k);   fp[2]:= pt(i+1,k+1);  fp[3]:= pt(i,k+1);
+               fe[1]:= e3(i,k);   fe[2]:= e1(i,k+1);    fe[3]:= e2(i,k);
+               end;  { with }
+          with edge[e1(i,k)] do
+               begin  ep1:= pt(i,k);  ep2:= pt(i+1,k); end;
+          with edge[e2(i,k)] do
+               begin  ep1:= pt(i,k);  ep2:= pt(i,k+1); end;
+          with edge[e3(i,k)] do
+               begin  ep1:= pt(i,k);  ep2:= pt(i+1,k+1); end;
+          end;  { for i }
+      end;  { for k }
+  for i:= 1 to n1-1 do
+      with edge[e1(i,n2)] do begin  ep1:= pt(i,n2); ep2:= pt(i+1,n2);  end;
+  for k:= 1 to n2-1 do
+      with edge[e2(n1,k)] do begin  ep1:= pt(n1,k); ep2:= pt(n1,k+1);  end;
+  if show_triangles then ne:= ne+nf;
+  for i:= 1 to ne do edge[i].vis:= true;
+  nf:= nf+nf;
+  for i:= 1 to nf do
+      with face[i] do
+           begin
+           plane_equ(p[fp[1]],p[fp[2]],p[fp[3]], nv,d,error);
+           npf:= 3; nef:= 3;
+           end;
+ end;  { aux_quadrangle_triang }
+{************}
+
+procedure aux_cylinder(n1,n2,np0,ne0,nf0: integer);
+var  n12,i,k : integer;   error : boolean;
+{***}
+function f(i,k :integer) : integer;  begin f := (k-1)*n1+i+nf0; end;
+function pt(i,k:integer) : integer;  begin pt:= (k-1)*n1+i+np0; end;
+function e1(i,k:integer) : integer;  begin e1:= (k-1)*n1+i+ne0; end;
+function e2(i,k:integer) : integer;  begin e2:= n12+(i-1)*(n2-1)+k+ne0; end;
+{****}
+begin
+  n12:= n1*n2;
+  for k:= 1 to n2 do
+      begin
+      for i:= 1 to n1 do
+          begin
+          with face[f(i,k)] do
+               begin
+               fp[1]:= pt(i,k);      fp[2]:= pt(i+1,k);
+               fp[3]:= pt(i+1,k+1);  fp[4]:= pt(i,k+1);
+               fe[1]:= e1(i,k);      fe[2]:= e2(i+1,k);
+               fe[3]:= e1(i,k+1);    fe[4]:= e2(i,k);
+               npf:= 4; nef:= 4;
+               end;  { with }
+          with edge[e1(i,k)] do
+               begin  ep1:= pt(i,k);  ep2:= pt(i+1,k);  vis:= true;  end;
+          if k<n2 then
+          with edge[e2(i,k)] do
+               begin  ep1:= pt(i,k);  ep2:= pt(i,k+1);  vis:= true;  end;
+          end;  { for i }
+      end;  { for k }
+  for k:= 1 to n2 do    { Korrektur }
+      begin
+      with face[f(n1,k)] do
+           begin
+           fp[2]:= pt(1,k);  fp[3]:= pt(1,k+1);
+           fe[2]:= e2(1,k);
+           end; { with }
+      edge[e1(n1,k)].ep2:= pt(1,k);
+      end;  { for k }
+  np:= np0+n12; ne:= ne0+2*n12-n1;  nf:= nf0+n12-n1;
+  for i:= nf0+1 to nf do
+      with face[i] do plane_equ(p[fp[1]],p[fp[2]],p[fp[3]], nv,d,error);
+end;  { aux_cylinder }
+{************}
+
+procedure aux_torus(n1,n2,np0,ne0,nf0: integer);
+ var  n12,i,k : integer;
+{***}
+function f(i,k :integer) : integer;  begin f := (k-1)*n1+i+nf0; end;
+function pt(i,k:integer) : integer;  begin pt:= (k-1)*n1+i+np0; end;
+function e1(i,k:integer) : integer;  begin e1:= (k-1)*n1+i+ne0; end;
+function e2(i,k:integer) : integer;  begin e2:= n12+(i-1)*n2+k+ne0; end;
+{****}
+ begin
+   n12:= n1*n2;
+   for k:= 1 to n2 do
+       begin
+       for i:= 1 to n1 do
+           begin
+           with face[f(i,k)] do
+                begin
+                fp[1]:= pt(i,k);      fp[2]:= pt(i+1,k);
+                fp[3]:= pt(i+1,k+1);  fp[4]:= pt(i,k+1);
+                fe[1]:= e1(i,k);      fe[2]:= e2(i+1,k);
+                fe[3]:= e1(i,k+1);    fe[4]:= e2(i,k);
+                npf:= 4; nef:= 4;
+                end;  { with }
+           with edge[e1(i,k)] do
+                begin  ep1:= pt(i,k);  ep2:= pt(i+1,k);  vis:= true;  end;
+           with edge[e2(i,k)] do
+                begin  ep1:= pt(i,k);  ep2:= pt(i,k+1);  vis:= true;  end;
+           end;  { for i }
+       end;  { for k }
+{ Korrektur: }
+   for k:= 1 to n2 do
+       begin
+       with face[f(n1,k)] do
+            begin
+            fp[2]:= pt(1,k);  fp[3]:= pt(1,k+1);
+            fe[2]:= e2(1,k);
+            end; { with }
+       with edge[e1(n1,k)] do ep2:= pt(1,k);
+       end;  { for k }
+   for i:= 1 to n1 do
+       begin
+       with face[f(i,n2)] do
+            begin
+            fp[3]:= np0+i+1;  fp[4]:= np0+i;
+            fe[3]:= ne0+i;
+            end;  { with }
+       with edge[e2(i,n2)] do ep2:= pt(i,1);
+       end;  { for i }
+   with face[f(n1,n2)] do
+        begin
+        fp[3]:= np0+1;  fp[4]:= np0+n1;    fe[3]:= ne0+n1;
+        end;
+  np:= np0+n12;  ne:= ne0+n12+n12;   nf:= nf0+n12;
+  for i:= nf0+1 to nf do
+      with face[i] do plane_equ(p[fp[1]],p[fp[2]],p[fp[3]], nv,d,error);
+ end;  { aux_torus }
+{************}
+
+procedure is_line_convex_polygon(p1,p2 : vt2d; p_pol : vts2d_pol; np : integer;
+                                           var t1,t2 : real; var ind : integer);
+{Berechnet die Parameter t1,t2 der Schnittpunkte der Strecke p1,p2
+ mit dem konvexen Polygon p_pol[0],...p_pol[np].
+ ind=0 bzw. 2 : Strecke  innerhalb bzw. ausserhalb, ind=1: sonst.
+ vts2d_pol: array[0..npfmax] of vt2d. }
+ var x1,y1,x2,y2,x21,y21,s,t,det,xi,yi,xi1,yi1 : real;  i,ns : integer;
+ begin
+   x1:= p1.x;  y1:= p1.y;  x2:= p2.x;  y2:= p2.y;
+   x21:= x2-x1;   y21:= y2-y1;    ns:= 0;
+   t1 :=0;  t2:=0;  i:=0;         xi:= p_pol[0].x;  yi:= p_pol[0].y;
+   while (ns<2) and (i<np) do
+         begin
+         xi1:= p_pol[i+1].x;  yi1:= p_pol[i+1].y;
+         det:= (xi1-xi)*y21 - x21*(yi1-yi);
+         if abs(det)>eps6 then
+            begin
+            s:= (x21*(yi-y1) - y21*(xi-x1))/det;
+            if (-eps7<=s) and (s<=1+eps7) then
+               begin
+               t:= ((yi-y1)*(xi1-xi)-(xi-x1)*(yi1-yi))/det;
+               if ns=0 then
+                  begin
+                  t1:= t;      ns:= 1;
+                  end
+                 else  {ns=1}
+                  begin
+                  t2:= t;   if abs(t2-t1)>eps6 then ns:= 2;
+                  end;
+               end;  { if }
+            end;  { if }
+         xi:= xi1;  yi:= yi1;  i:= i+1;
+         end; { while }
+   ind:= 2;
+   if ns=2 then
+      begin
+      ind:= 1;   {Strecke schneidet wenigstens eine Seite}
+      if t2<t1 then  change1d(t1,t2);
+      if (t1<=0) and (t2>=1) then ind:= 0;  {Strecke innerh.}
+      if (t1>=1) or  (t2<=0) then ind:= 2;  {Strecke ausserh.}
+      end;
+ end;  { is_line_convex_polygon }
+{*************}
+
+procedure intmint(a,b,c,d: real; var e1,f1,e2,f2: real; var ind: integer);
+{Berechnet die Intervall-Differenz  [a,b] \ [c,d] .
+ ind=0: leer, ind=1: 1 Interv., ind=2: 2 Interv.}
+ var aa,bb,cc,dd : real;
+ begin
+   aa:= min(a,b);   bb:= max(a,b);
+   cc:= min(c,d);   dd:= max(c,d);
+   e1:= aa;         f1:= bb;        ind:= 1;
+   if (cc<=aa) and (dd>=bb)   then  ind:=0
+     else
+      begin
+      if (cc<=aa) and (dd>aa)  then  e1:= dd;
+      if (dd>=bb) and (cc<bb)  then  f1:= cc;
+      if (cc>aa)  and (dd<bb)  then
+         begin
+         f1:= cc;  e2:= dd;  f2:= bb;  ind:= 2;
+         end;
+      end;
+ end; { intmint }
+{*************}
+
+procedure cp_lines_before_convex_faces(oriented_faces,is_permitted,newstyles : boolean);
+{Projiziert und zeichnet Kantenteile VOR (orientierten) ebenen n-Ecken.
+ oriented_faces=true: die Flaechen sind orientiert,
+ is_permitted=true: Kanten duerfen die Flaechen schneiden.
+ Aus dem Hauptprogramm muessen bereitstehen:
+ np, ne, nf: Anzahl der Punkte, Kanten, Flaechen,
+ p[1],...,p[np] : Punkte,
+ face[i].fp[k] (face[i].fe[k]): k-ter Punkt (k-te Kante) in i-ter Flaeche,
+ face[i].npf (face[i].nef): Anzahl der Punkte (Kanten) in der i-ten Flaeche,
+ edge[i].ep1 (edge[i].ep2): Anfangs-(End-)Punkt der i-ten Kante,
+ face[i].nv,face[i].d: Koeffizienten der Ebenengleichung.}
+ label 5,10;
+ var  t1,t2,tt,dispt,xemin,xemax,yemin,yemax,zemin,d1,d2,dd1,dd2,a,b,
+      tt1,tt2,ts,disp1,disp2,dispt1,dispt2,test1,test2 : real;
+      i,j,k,l,m,i1,i2,nseg,nseg0,ind : integer;
+      p1,p2,pd,pt, pt1,pt2,ps : vt3d;
+      pp1,pp2,ppd,qq1,qq2,pps : vt2d;
+      pp : vts2d; 
+      par1,par2,pa1,pa2,pa3,pa4 : r_array_seg;  {r_array_seg : s. HP    }
+      p_pol : vts2d_pol;                        {vts2d_pol   : s. HP    }
+      nt : i_array_seg;                         {i_array_seg : s. HP    }
+{******}
+ begin
+    if nf>100 then
+    begin
+       writeln; writeln('####:  wait !!! (for hiddenline-alg.)  '); writeln;
+    end;
+    if oriented_faces then
+      begin
+      for i:= 1 to nf do face[i].vis:= false;
+      for i:= 1 to ne do edge[i].vis:= false;
+      end;
+{ Koordinaten der Bildpunkte: }
+   cp_vts3d_vts2d_spez(p,1,np, pp,pdist);
+{ Fenster und Normalentest fuer die Flaechen: }
+   for i:= 1 to nf do
+       begin
+       with face[i] do
+            begin
+            discentre:= scalarp3d(nv,centre) - d;
+            if ((discentre>=0) or (not oriented_faces)) then
+               begin
+               if oriented_faces then
+                  begin      { Normalentest: Flaeche ist sichtbar }
+                  vis:= true;
+                  for k:= 1 to nef do edge[fe[k]].vis:= true;
+                  end;
+               with box do
+                    begin
+                    with pp[fp[1]] do   { Anfangswerte }
+                         begin xmin:= x; xmax:= x;  ymin:= y;  ymax:= y; end;
+                    zmax:= pdist[fp[1]];
+                    for k:= 2 to npf do
+                        begin
+                        with pp[fp[k]] do     { Flaechenfenster }
+                             begin
+                             xmin:= min(xmin,x);  ymin:= min(ymin,y);
+                             xmax:= max(xmax,x);  ymax:= max(ymax,y);
+                             end;  { with }
+                        zmax:= max(zmax,pdist[fp[k]]);
+                        end;  { for k }
+                    end;  { with box }
+               end;  { if }
+            end; { with  face }
+       end; { for i }
+{ Test und Zeichnen der Kanten(-teile):}
+   for i:= 1 to ne do                           { Beginn der KANTENschleife }
+       begin
+       if not edge[i].vis then goto 10;
+       par1[1]:= 0;  par2[1]:= 1;  nseg:= 1;    { : 1 sichtb. Anfangsintervall}
+  { Punkte und Fenster der Kante : }
+       i1:= edge[i].ep1;  i2:= edge[i].ep2;
+       p1:= p[i1];  p2:= p[i2];  pp1:= pp[i1]; pp2:= pp[i2];
+       xemax:= max(pp1.x,pp2.x);   yemax:= max(pp1.y,pp2.y);
+       xemin:= min(pp1.x,pp2.x);   yemin:= min(pp1.y,pp2.y);
+       zemin:= min(pdist[i1],pdist[i2]);
+       for j:= 1 to nf do                       { Beginn der FLAECHENschleife }
+           with face[j] do
+           begin
+           if not vis and oriented_faces then goto 5;
+  { Fenstertest mit j-ter Flaeche: }
+           if (xemax<=box.xmin) or (xemin>=box.xmax)
+             or (yemax<=box.ymin) or (yemin>=box.ymax) or (box.zmax<=zemin)
+               then   goto 5; { Kante wird nicht von j-ter Flaeche verdeckt }
+  { Kante i Kante von j-ter Flaeche ?:}
+           for k:= 1 to nef do  if i=fe[k] then goto 5; { naechste Flaeche}
+  { Schnitt der Kante mit dem (konvexen) Flaechenpolygon in der Bildtafel: }
+           for k:= 1 to npf do p_pol[k]:= pp[fp[k]];  p_pol[0]:= p_pol[npf];
+           is_line_convex_polygon(pp1,pp2,p_pol,npf,t1,t2,ind);
+           if ind=2 then  goto 5 ;   {Kante nicht verdeckt=> naechste Flaeche }
+	   d1:= distance3d(p1,centre);     d2:= distance3d(p2,centre);
+	   with pp1 do dd1:=sqrt(x*x+y*y+sqr(distance)); {fuer Testpunkte}
+	   with pp2 do dd2:=sqrt(x*x+y*y+sqr(distance)); 
+	   a:= d1/dd1;  b:= d2/dd2;
+           if not is_permitted then  {Kante darf die Flaeche NICHT durchdr.}
+  { Testpunkt und Test, ob Kante vor oder hinter der Flaeche : }
+              begin
+              tt:= ( max(0,t1) + min(1,t2) )/2;     diff3d(p2,p1, pd);
+	      tt:= a*tt/(b+tt*(a-b)); {Korrekt. (Zentralproj. keine lin. Abb.)}
+              lcomb2vt3d(1,p1, tt,pd, pt);
+              dispt:= scalarp3d(nv,pt) - d;    {Distanz Testpunkt-Ebene}
+              if (dispt*discentre>0)  then   goto 5;
+                { Kante vor der Ebene => naechste Flaeche }
+              end;  { if not is_permitted }
+           if is_permitted then  { Kante darf die Flaeche durchdringen }
+ { Bestimmung des Teils der Kante, der hinter der Flaeche liegt: }
+              begin
+              tt1:= max(0,t1);  tt2:= min(1,t2);   diff3d(p2,p1, pd);
+	      tt1:= a*tt1/(b+tt1*(a-b));{Korrekt. (Zentralproj. keine lin. Abb.)}
+	      tt2:= a*tt2/(b+tt2*(a-b));{                "                      }
+              lcomb2vt3d(1,p1,tt1,pd, pt1);      lcomb2vt3d(1,p1,tt2,pd, pt2);
+              dispt1:= scalarp3d(nv,pt1) - d;    {Distanz 1.Testpunkt-Ebene}
+              dispt2:= scalarp3d(nv,pt2) - d;    {Distanz 2.Testpunkt-Ebene}
+              test1:= dispt1*discentre;          test2:= dispt2*discentre;
+              if ((test1>=0) and (test2>=0))  then   goto 5;
+             { Kante vor der Ebene => naechste Flaeche }
+              if dispt1*dispt2<0 then
+                 begin
+                 disp1:= scalarp3d(nv,p1) - d;   disp2:= scalarp3d(nv,p2) - d;
+                 ts:= disp1/(disp1-disp2);
+                 lcomb2vt3d(1,p1,ts,pd, ps);  cp_vt3d_vt2d(ps,pps);
+                 diff2d(pp2,pp1, ppd); {Zentralproj. ist keine lin. Abb. !!!}
+                 ts:= (scalarp2d(pps,ppd)-scalarp2d(pp1,ppd))/scalarp2d(ppd,ppd);
+                 if test1<0 then t2:= ts else t1:= ts;
+                 end;
+              end; { if is_permitted}
+              for l:= 1 to nseg do   {weiterhin sichtbare Intervalle:}
+                  intmint(par1[l],par2[l],t1,t2,
+                                         pa1[l],pa2[l],pa3[l],pa4[l],nt[l]);
+              nseg0:= nseg; m:= 0;
+              for l:= 1 to nseg0 do          { neue Intervallteilung }
+                  begin
+                  if nt[l]=0 then  nseg:= nseg-1   { 1 Segment weniger }
+                    else
+                     begin
+                     m:= m+1;
+                     par1[m]:= pa1[l];  par2[m]:= pa2[l];
+                     if nt[l]=2 then      { 1 Segment mehr }
+                        begin
+                        m:= m+1;  nseg:= nseg+1;
+                        par1[m]:= pa3[l];  par2[m]:= pa4[l];
+                        end;
+                     end;  { if }
+                  end; { for l }
+          if nseg<1 then goto 10;
+5:        end; { with }
+       for k:=1 to nseg do   { Zeichnen der sichtb. Segmente der i-ten Kante}
+           begin
+           diff2d(pp2,pp1, ppd);
+           lcomb2vt2d(1,pp1, par1[k],ppd, qq1);
+           lcomb2vt2d(1,pp1, par2[k],ppd, qq2);
+	   if newstyles then
+	      with edge[i] do begin new_color(color); new_linewidth(linew); end;
+           line2d(qq1,qq2,0);
+           end;
+10:    end; { for i (Kantenschleife) }
+ end; { cp_lines_before_convex_faces }
+{*************}
+
+
+
+procedure is_interv_interv(var a,b,c,d,aa,bb : real; var inters: boolean);
+{Berechnet den Schnitt der Intervalle [a,b], [c,d] .}
+var a1,a2,b1,b2 : real;
+begin
+  a1:= min(a,b);    b1:= max(a,b);
+  a2:= min(c,d);    b2:= max(c,d);
+  aa:= max(a1,a2);  bb:= min(b1,b2);
+  if bb<=aa then inters:= false else inters:= true;
+end; {is_interv_interv}
+{*****}
+
+procedure box3d_of_pts(var p : vts3d_pol; np: integer; var box : box3d_dat);
+var   i: integer;
+{Bestimmt den zugehoerigen (kleinsten) achsenparallelen Quader.}
+begin
+  with box do
+       begin
+       with p[1] do
+            begin
+            xmin:= x; xmax:= x;   ymin:= y; ymax:= y;  zmin:= z; zmax:= z;
+            end;
+       for i:= 2 to np do with p[i] do
+           begin
+           xmin:= min(xmin,x);  xmax:= max(xmax,x);
+           ymin:= min(ymin,y);  ymax:= max(ymax,y);
+           zmin:= min(zmin,z);  zmax:= max(zmax,z);
+           end;  { for }
+      end;  { with }
+end;  { box3d_of_pts }
+{************}
+(*
+function is_two_boxes3d(var box1,box2 : box3d_dat) : boolean;
+{Schnitt zweier achsenparalleler Quader.}
+begin
+  is_two_boxes3d:= false;
+  if (box1.xmax>box2.xmin) then if (box1.xmin<box2.xmax) then
+     if (box1.ymax>box2.ymin) then if (box1.ymin<box2.ymax) then
+        if (box1.zmax>box2.zmin) then if (box1.zmin<box2.zmax) then
+           is_two_boxes3d:= true;
+end; {is_two_boxes3d }
+*)
+{*************}
+
+
+function is_two_boxes3d(var box1,box2 : box3d_dat) : boolean;
+{Schnitt zweier achsenparalleler Quader.}
+begin
+
+  if ((box1.xmax>box2.xmin) and (box1.xmin<box2.xmax) and
+      (box1.ymax>box2.ymin) and (box1.ymin<box2.ymax) and
+      (box1.zmax>box2.zmin) and (box1.zmin<box2.zmax))
+     then
+       is_two_boxes3d:=true
+     else
+       is_two_boxes3d:=false;
+end;  {is_two_boxes3d }
+{*************}
+
+procedure is_line_conv_pol_in_plane3d(var pl,rl: vt3d; var pp : vts3d_pol;
+                                      npp : integer;
+                                      var t1,t2 : real; var inters : boolean);
+{Schnitt eines konv. Polygons mit einer in der Polygonebene liegenden Gerade.}
+var pp1,v1,v2,pl2 : vt3d;      qql1,qql2 : vt2d;
+    qq : vts2d_pol;    error : boolean;   i,ind : integer;
+begin
+  pp1:= pp[1];  diff3d(pp[npp],pp1, v2);  diff3d(pp[2],pp1, v1);
+  for i:= 3 to npp-1 do
+      ptco_plane3d(pp1,v1,v2,pp[i], qq[i].x,qq[i].y,error);
+  put2d(0,0, qq[1]);    put2d(1,0, qq[2]);
+  put2d(0,1, qq[npp]);  qq[0]:= qq[npp];
+  ptco_plane3d(pp1,v1,v2,pl,qql1.x,qql1.y,error);
+  sum3d(pl,rl, pl2);
+  ptco_plane3d(pp1,v1,v2,pl2, qql2.x,qql2.y,error);
+  is_line_convex_polygon(qql1,qql2,qq,npp, t1,t2,ind);
+  if t1<>t2 then inters:= true else inters:= false;
+end;  { is_line_conv_pol_in_plane3d }
+{***********}
+procedure is_n1gon_n2gon3d(var pp1,pp2: vts3d_pol; np1,np2: integer;
+                           var  ps1,ps2 : vt3d; var intersection : boolean);
+{Berechnet die Schnittstrecke zweier ebener konvexer Polygone im Raum, die
+ NICHT in einer Ebene liegen.}
+var   box1,box2 : box3d_dat;    t1,t2,s1,s2,s,t : real;
+      nv1,nv2,ps,rs : vt3d;  d1,d2 : real;
+      inters1,inters2,error,inters : boolean;
+begin
+  intersection:= false;
+  box3d_of_pts(pp1,np1, box1);  box3d_of_pts(pp2,np2, box2);
+  if is_two_boxes3d(box1,box2) then
+     begin
+     plane_equ(pp1[1],pp1[2],pp1[3], nv1,d1,error);
+     plane_equ(pp2[1],pp2[2],pp2[3], nv2,d2,error);
+     is_plane_plane(nv1,d1,nv2,d2, ps,rs, error);
+     if not error then
+        begin
+        is_line_conv_pol_in_plane3d(ps,rs,pp1,np1, t1,t2,inters1);
+        is_line_conv_pol_in_plane3d(ps,rs,pp2,np2, s1,s2,inters2);
+        if inters1 and inters2 then
+           begin
+           is_interv_interv(t1,t2,s1,s2, s,t,inters);
+           if inters then
+              begin
+              lcomb2vt3d(1,ps,s,rs, ps1);
+              lcomb2vt3d(1,ps,t,rs, ps2);
+              intersection:= true;
+              end;   { if inters }
+           end;  { if inters1,inters2 }
+        end;  { if not error }
+     end; { if is_boxes }
+end;  { is_n1gon_n2gon }
+{**************}
+
+procedure boxes_of_faces;
+var i,j : integer;  pp : vts3d_pol;
+begin
+  for i:= 1 to nf do
+      with face[i] do
+           begin
+           for j:= 1 to npf do pp[j]:= p[fp[j]];
+           box3d_of_pts(pp,npf, box);
+           end;
+end;  { boxes_of_faces }
+{************}
+
+procedure is_face_face(i,k: integer;  var ps1,ps2 : vt3d;
+                                      var intersection: boolean);
+{Berechnet die Schnittstrecke zweier nicht in einer Ebene liegenden
+ (konvexen) Flaechen eines Polyeders.}
+var  t1,t2,s1,s2,s,t : real;      ps,rs : vt3d;  ppf: vts3d_pol;
+     error,inters1,inters2,inters : boolean;  j: integer;
+begin
+  intersection:= false;
+  if is_two_boxes3d(face[i].box,face[k].box) then
+     begin
+     is_plane_plane(face[i].nv,face[i].d,face[k].nv,face[k].d, ps,rs, error);
+     if not error then
+        begin
+        with face[i] do for j:= 1 to npf do ppf[j]:= p[fp[j]];
+        is_line_conv_pol_in_plane3d(ps,rs,ppf,face[i].npf, t1,t2,inters1);
+        with face[k] do for j:= 1 to npf do ppf[j]:= p[fp[j]];
+        is_line_conv_pol_in_plane3d(ps,rs,ppf,face[k].npf, s1,s2,inters2);
+        if inters1 and inters2 then
+           begin
+           is_interv_interv(t1,t2,s1,s2, s,t,inters);
+           if inters then
+              begin
+              lcomb2vt3d(1,ps,s,rs, ps1);
+              lcomb2vt3d(1,ps,t,rs, ps2);
+              intersection:= true;
+              end;   { if inters }
+           end;  { if inters1,inters2 }
+        end;  { if not error }
+     end; { if is_boxes }
+end;  { is_face_face }
+{**************}

cdg0gv-2010/tools/Makefile

@@ -0,0 +1,22 @@
+#PC=	ppc386
+PC=	ppcx64
+
+CDGDIR=        ..
+INCLDIR=        $(CDGDIR)/include
+UNITDIR=        $(CDGDIR)/units
+TOOLDIR=        $(CDGDIR)/tools
+
+PFLAGS=		-So -Fu../units -I../include
+
+%:		%.p 
+		$(PC) $(PFLAGS)  $<
+
+tools:		
+		$(PC) $(PFLAGS) pldv.p
+		rm -f *~ *.o 
+
+clean:
+		rm -f *~ *.o pldv
+
+
+

cdg0gv-2010/tools/include

@@ -0,0 +1 @@
+../include

cdg0gv-2010/tools/pldv.p

@@ -0,0 +1,319 @@
+program pldv;
+uses  Unix {linux};
+const
+{$i include/geoconst.pas}
+
+type
+{$i include/geotype.pas}
+
+var
+{$i include/geovar.pas}
+
+const
+  AutoScale    : Boolean = false;
+  Frame        : Boolean = false;
+  Fixed_X_Size : Boolean = false;
+  Fixed_Y_Size : Boolean = false;
+  Postscript   : Boolean = false;
+  eps          : Boolean = false;
+  mono	       : Boolean = false;
+  RandDefault  : Boolean = true;  
+  
+var
+  breite,hoehe,
+  xnullp,ynullp,
+  xmin,xmax,
+  ymin,ymax,r,
+  x_size,y_size,
+  faktor,
+  Rand,
+  x1,y1,x2,y2        : Real;
+  anz,
+  ind,i,idummy       : Integer;
+  C                  : Char;
+  PLDFileName,
+  PSFileName         : string;
+  PLDFILE            : text;
+  s                  : String;
+  color              : integer;
+  Option_Error,
+  anfang,newname     : Boolean;
+  p                  : vts2d;
+
+
+{$i include/geoproc.pas}
+{$i include/proc_ag.pas}
+{i postscr.pas}
+
+
+procedure remapx(var x:real);
+begin
+  x := (x/pt_per_mm + xnullp) * faktor + Rand;
+end;
+
+procedure remapy(var y:real);
+begin
+  y := (y/pt_per_mm + ynullp) * faktor + Rand;
+end;
+
+{******************}
+
+begin   {main}
+
+  PLDFilename:='';
+  PSFilename:='geodummy.eps';
+  faktor:=1.0;  
+{  Postscript:= true;}
+  eps       := true;
+  newname   := false;
+
+  { Parser for options }
+
+  for i := 1 to ParamCount do
+  begin
+    Option_Error := False;
+    s := ParamStr(i);
+    if length(s) > 0 then
+    begin
+      Case s[1] of
+         '-' : begin
+                if length(s) > 1 then
+                Case s[2] of
+                  'o' : begin
+                          s := Copy(s,3,length(s)-2);
+        		  PSFileName:=s; newname:= true;
+                        end;
+                  'p' : begin
+		          Postscript:= true;
+                          s := Copy(s,3,length(s)-2);
+ 		          if s='ps' then eps:= false
+                        end;
+                  'z' : begin
+                          s := Copy(s,3,length(s)-2);
+                          Val(s,idummy,ind);
+                          if ind <> 0 then begin
+                            faktor := 1;
+                            Option_Error := True;
+                          end else faktor := idummy / 100.0;
+                        end;
+                  'a' : AutoScale := True;
+                  'm' : mono := True;
+                  'x' : begin
+                          Fixed_X_Size := True;
+                          Fixed_Y_Size := false;
+                          s := Copy(s,3,length(s)-2);
+                          Val(s,idummy,ind);
+                          if ind <> 0 then begin
+                            Fixed_X_size := false;
+                            Option_Error := True;
+                          end else x_size := idummy *1.0;
+                        end;
+                  'y' : begin
+                          Fixed_Y_Size := True;
+                          Fixed_X_Size := false;
+                          s := Copy(s,3,length(s)-2);
+                          Val(s,idummy,ind);
+                          if ind <> 0 then begin
+                            Fixed_Y_size := false;
+                            Option_Error := True;
+                          end else y_size := idummy *1.0;
+                        end;
+                  'r' : begin
+                          RandDefault := False;
+                          s := Copy(s,3,length(s)-2);
+                          Val(s,idummy,ind);
+                          if ind <> 0 then begin
+                            RandDefault := True;
+                            Option_Error := True;
+                          end else Rand := idummy *1.0;
+                        end;
+                  else Option_Error := True; { nur -  }
+                 end;
+              end; 
+          else PLDFileName:=s;
+      end { Case }
+    end; {if}
+
+    if Option_Error then
+    begin
+      Writeln('inadmissible option : ', ParamStr(i) );
+      halt;
+    end;
+  end;  { For }
+
+  if PLDFileName='' then begin
+    writeln;
+    writeln('PLDV by Andreas Goerg / Erich Hartmann (Version: 10.12.2002)');
+    writeln('for PLD-files version >=1.0 ');
+    writeln('call: PLDV {PLD-file-name} [-p{ps,eps}] [-o{file}] [-z{%}] [-a]');
+    writeln('                                 [-x{width}] [-y{height}] [-r{border}] ');
+    writeln(' -p : a ps/eps-file is generated with the same name as the pld-file');
+    writeln('      or the name specified in -o.... ');
+    writeln(' -z : zoom (per cent) ');
+    writeln(' -a : adjust the drawing area ');
+    writeln(' -x : drawing scaled to {width} mm ');
+    writeln(' -y : drawing scaled to {height} mm ');
+    writeln(' -m : (monochrome) ignoring colors');
+{    writeln(' -q : keine Grafikausgabe');}      
+    writeln(' -r : border width= {border} mm around the drawing'); 
+    halt;
+  end;  
+
+   if (not newname) and (Postscript) then
+   begin
+       s:= Copy(pldfilename,1,length(pldfilename)-3);
+       if eps then s:= s+'eps' else s:= s+'ps';
+       psfilename:= s; 
+   end; 
+
+  writeln('PLDV: A Postscript-file ',PSFileName, ' is generated');
+    
+  assign(PLDFile,PLDFilename);
+  reset(PLDFile);
+
+  Readln(PLDFile,s);
+  if s <> 'Ver1.0' then 
+  begin
+    Writeln('This PLD-File is not of version >= Ver 1.0');
+    halt;
+  end;  
+  
+  xnullp := 0;
+  ynullp := 0;
+
+  if AutoScale or Fixed_X_Size or Fixed_Y_Size then
+  begin
+    Read(PLDFile,breite,hoehe);
+    anfang := True;
+    repeat
+      Read(PLDFile,C);
+      Case C of
+      'P' : begin
+              Read(PLDFile,x1,y1,ind);
+              if anfang then
+              begin
+                xmax := x1; xmin := xmax;
+                ymax := y1; ymin := ymax;
+                anfang := false;
+              end
+              else
+              begin
+                if x1 > xmax then xmax := x1;
+                if x1 < xmin then xmin := x1;
+                if y1 > ymax then ymax := y1;
+                if y1 < ymin then ymin := y1;
+              end;
+            end;
+      'L' : begin
+              Read(PLDFile,x1,y1,x2,y2,ind);
+              if anfang then
+              begin
+                xmax := x1; xmin := xmax;
+                ymax := y1; ymin := ymax;
+                anfang := false;
+              end
+              else
+              begin
+                if x1 > xmax then xmax := x1;
+                if x1 < xmin then xmin := x1;
+                if y1 > ymax then ymax := y1;
+                if y1 < ymin then ymin := y1;
+                if x2 > xmax then xmax := x2;
+                if x2 < xmin then xmin := x2;
+                if y2 > ymax then ymax := y2;
+                if y2 < ymin then ymin := y2;
+              end;
+            end;
+      'K' : begin
+              Read(PLDFile,anz,ind);
+              for i:= 1 to anz do 
+              begin
+                Read(PLDFile,x1,y1);
+                if anfang then
+                begin
+                  xmax := x1; xmin := xmax;
+                  ymax := y1; ymin := ymax;
+                  anfang := false;
+                end
+                else
+                begin
+                  if x1 > xmax then xmax := x1;
+                  if x1 < xmin then xmin := x1;
+                  if y1 > ymax then ymax := y1;
+                  if y1 < ymin then ymin := y1;
+                end;
+              end;                      
+            end;
+    end;
+    until (C = '*') or EOF(PLDFile);
+
+    Reset(PLDFile);
+    Readln(PLDFile,s);
+
+    Breite := (xmax-xmin)/pt_per_mm;
+    Hoehe  := (ymax-ymin)/pt_per_mm;
+    xnullp := -xmin/pt_per_mm;
+    ynullp := -ymin/pt_per_mm;
+
+    if RandDefault then Rand := 3;
+
+  end
+  else
+  begin
+    if RandDefault then Rand := 0;
+  end; 
+
+  if Fixed_Y_Size then  faktor := pt_per_mm*y_size/(ymax-ymin);
+  if Fixed_X_Size then  faktor := pt_per_mm*x_size/(xmax-xmin);
+
+   ps_graph_on(PSFilename,eps);
+  
+  if AutoScale or Fixed_X_size or Fixed_Y_Size then Read(PLDFile,r,r)
+                                               else Read(PLDFile,breite,hoehe);
+    ps_Draw_Area(breite*faktor+2*Rand,hoehe*faktor+2*Rand,0,0,eps);  
+
+  repeat
+    Read(PLDFile,C);
+    Case C of
+    'P' : begin
+            Read(PLDFile,x1,y1,ind); remapx(x1); remapy(y1);
+            ps_pointc2d(x1,y1,ind);
+          end;
+    'L' : begin
+            Read(PLDFile,x1,y1,x2,y2,ind); remapx(x1); remapy(y1); remapx(x2); remapy(y2);
+            ps_linec2d(x1,y1,x2,y2,ind);
+          end;
+    'K' : begin
+            Read(PLDFile,anz,ind);
+            for i:= 1 to anz do 
+            begin
+              Read(PLDFile,x1,y1); remapx(x1); remapy(y1);
+              put2d(x1,y1,p[i-1]);
+            end;
+            ps_curve2d(p,0,anz-1,ind);
+            end;
+    'C' : begin
+            Read(PLDFile,color);
+            if not mono then
+            begin
+              ps_new_color(color);
+            end;
+          end;
+    'W' : begin
+            Read(PLDFile,i);
+            ps_new_linewidth(i/10.0);
+          end;
+    end;
+  until (C = '*') or EOF(PLDFile);
+
+  ps_Draw_end;            
+  ps_Graph_Off;  
+
+{Show psfile:}
+   s:= 'gv '; s:= s+psfilename;  
+   writeln('### display: ',s);
+   {* shell(s); *}
+   fpSystem(s);
+
+end.
+

cdg0gv-2010/units/Makefile

@@ -0,0 +1,20 @@
+#PC=	ppc386
+PC=	ppcx64
+
+CDGDIR=        ..
+INCLDIR=        $(CDGDIR)/include
+UNITDIR=        $(CDGDIR)/units
+TOOLDIR=        $(CDGDIR)/tools
+
+PFLAGS=		-So -Fu../units -I../include
+
+%:		%.p 
+		$(PC) $(PFLAGS)  $<
+units:		
+		$(PC) $(PFLAGS)  geograph.p 
+		$(PC) $(PFLAGS)  hiddenl.p 
+clean:
+		rm -f *~ *.o *.ppu 
+
+
+

cdg0gv-2010/units/geograph.p

@@ -0,0 +1,36 @@
+UNIT geograph;
+
+{SMARTLINK ON}
+
+INTERFACE
+
+USES Unix {linux};
+{************************}
+{***  g e o g r a p h ***}
+{************************}
+
+const {$i include/geoconst.pas}
+type  {$i include/geotype.pas}
+var   {$i include/geovar.pas}
+
+{$i include/head_geo.pas}
+{$i include/head_ag.pas}
+{$i include/head_pp.pas}
+{ i include/head_ks.pas}
+{ i include/head_pks.pas}
+{ i include/head_qua.pas}
+{ i include/head_pqu.pas}
+{$i include/head_zp.pas}
+
+IMPLEMENTATION
+
+{$i include/geoproc.pas}
+{$i include/proc_ag.pas}
+{$i include/proc_pp.pas}
+{ i include/proc_ks.pas}
+{ i include/proc_pks.pas}
+{ i include/proc_qua.pas}
+{ i include/proc_pqu.pas}
+{$i include/proc_zp.pas}
+
+END. {of IMPLEMENTATION}

cdg0gv-2010/units/hiddenl.p

@@ -0,0 +1,47 @@
+UNIT hiddenl;
+
+INTERFACE
+
+USES GEOGRAPH;
+{************************}
+{***  h i d d e n l ***}
+{************************}
+
+{fuer Hiddenline:}
+
+const {Achtung: es muss  array_size>=nfmax  sein  !!!}
+       nfmax= 10000;  nemax=20000;    nsegmax=10;  npfmax=10;
+
+type   vts2d_pol  = array[0..npfmax] of vt2d;
+       vts3d_pol  = array[0..npfmax] of vt3d;
+       r_array_seg  = array[0..nsegmax] of real;
+       i_array_seg  = array[0..nsegmax] of integer;
+       box3d_dat = record
+                   xmin,xmax,ymin,ymax,zmin,zmax : real;
+                   end;
+       face_dat = record
+                  npf,nef : integer;
+                  fp,fe : array[1..npfmax] of integer;
+                  vis : boolean;
+                  box : box3d_dat;
+                  discentre,d : real;
+                  nv : vt3d;
+                  end;
+       edge_dat = record
+                  vis : boolean;
+                  ep1,ep2,color,linew : integer;
+                  end;
+
+var    ne,nf,np: integer;  {Anzahl der Kanten, Facetten,Punkte}
+       p    : vts3d;       {Punkte des Polyeders} 
+       face : array[1..nfmax] of face_dat;
+       edge : array[1..nemax] of edge_dat;
+       pdist: r_array;     {pdist[i]: Abstand d. i-ten Punktes von d. Bildeb.} 
+       error,oriented_faces,is_permitted,newstyles: boolean;
+{$i include/head_zpo.pas}
+
+IMPLEMENTATION
+
+{$i include/proc_zpo.pas}
+
+END. {of IMPLEMENTATION}

cdg0gv-2010/units/include

@@ -0,0 +1 @@
+../include