summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/dvips/pst-3dplot/pst-3dplot.pro
blob: 781251deea2b98f74c31ceee3b04d77f02da66ae (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
%% $Id: pst-3dplot.pro 298 2010-03-13 08:46:53Z herbert $
%%
%% This is file `pst-3dplot.pro',
%%
%% IMPORTANT NOTICE:
%%
%% Package `pst-3dplot.tex'
%%
%% Herbert Voss <voss _at_ PSTricks.de>
%%
%% This program can be redistributed and/or modified under the terms
%% of the LaTeX Project Public License Distributed from CTAN archives
%% in directory macros/latex/base/lppl.txt.
%%
%% DESCRIPTION:
%%   `pst-3dplot' is a PSTricks package to draw 3d curves and graphical objects
%%
%%
%% version 0.31 / 2010-02-20  Herbert Voss <hvoss _at_ tug.org>
%% with contributions of Darrell Lamm <darrell.lamm _at_ gtri.gatech.edu<
%%            
%
/tx@3DPlotDict 200 dict def
tx@3DPlotDict begin
%
/printDot { gsave 2 copy 2 0 360 arc fill stroke grestore } def
%
/saveCoor { 
  dzUnit mul /z ED
  dyUnit mul /y ED
  dxUnit mul /x ED
} def
%
/3Dto2D { % true or false on stack
  { RotatePoint } if
  1 { %  dummy loop, will run only 1 time, allows exit 
    coorType 0 le {                                               % the default |
      /x2D x leftHanded not { neg } if Alpha cos mul y Alpha sin mul add def %  /\  co system
      /y2D x leftHanded { neg } if Alpha sin mul y Alpha cos mul add neg Beta sin mul z Beta cos mul add def
      exit } if
    coorType 1 le { 
      /x2D y x Alpha 90 sub sin mul sub def  %  |/_  co system, no shortened x axis
      /y2D z x Alpha 90 sub cos mul sub def 
      exit } if
    coorType 2 le { % coorType |/_ with a 1/sqrt(2) shortend x-axis and 135 degrees 
      /x2D y x 0.5 mul sub def
      /y2D z x 0.5 mul sub def 
      exit } if
    coorType 3 le { % coorType |/_ with a 1/sqrt(2) shortend x-axis and 135 degrees 
      /x2D y x -0.5 mul sub def
      /y2D z x -0.5 mul sub def 
      exit } if
    coorType 4 le { % Normalbild in Trimetrie Skalierung so, dass coorType2
       /x2D x -0.5 mul y 1 mul add def
       /y2D x -0.5 mul y -0.25 mul add z 1 mul add def
       exit } if
  } repeat
} def
/ConvertTo2D { true 3Dto2D } def
/ConvertTo2DWithoutRotating { false 3Dto2D } def
%
/Conv3D2D { /z ED /y ED /x ED ConvertTo2D x2D y2D } def
%
/ConvertToCartesian {
  /latitude exch def
  /longitude exch def
  /Radius exch def
  1 { %  dummy loop, will run only 1 time, allows exit
    SphericalCoorType 0 le {                                               % the default |
     /z { Radius latitude sin mul } def
     /x { Radius longitude cos mul latitude cos mul } def
     /y { Radius longitude sin mul latitude cos mul } def
      exit } if
    SphericalCoorType 2 le {
     /z { Radius longitude cos mul } def
     /x { Radius longitude sin mul latitude cos mul} def
     /y { Radius longitude sin mul latitude sin mul } def
      exit } if
  } repeat
} def
%
/ConvCylToCartesian { % r phi h -> x y z
  3 1 roll			% h r phi
  /Phi ED
  /Radius ED			% h->z on stack
  Radius Phi cos mul exch 	% x z
  Radius Phi sin mul exch	% x y z
} def
%
/SphericalTo2D {
  x y z ConvertToCartesian ConvertTo2D
} def
%
/CylinderTo2D { %  r phi h
  x y z ConvCylToCartesian ConvertTo2D
} def
%
/convertStackTo2D {
  counttomark
  /n ED /n3 n 3 div cvi def
  n3 {
    n -3 roll
    SphericalCoor { ConvertToCartesian } { saveCoor } ifelse
    ConvertTo2D
    x2D xUnit y2D yUnit
    /n n 1 sub def
  } repeat
} def
%
% the angle in the parameter equation for an ellipse is not proportional to the real angle!
% phi=atan(b*tan(angle)/a)+floor(angle/180+0.5)*180
%
/getPhi { % on stack: vecA vecB angle 
  3 dict begin
  /angle exch def /vecB exch def /vecA exch def
  angle cvi 90 mod 0 eq { angle } { vecA angle tan mul vecB atan 
  angle 180 div .5 add floor 180 mul add } ifelse 
  end
} def
%
/RotSet (set ) def
%
/eulerRotation false def
% Matrix multiplication procedure
/matmul {

  /M@tMulDict 20 dict def
  M@tMulDict begin
  /m2 ED
  /m1 ED
  m1 dup length 2 sub 2 getinterval aload pop
  /col1max ED
  /row1max ED
  m2 dup length 2 sub 2 getinterval aload pop
  /col2max ED
  /row2max ED
  /m3 row1max col2max mul 2 add array def
  m3 dup length 2 sub row1max col2max 2 array astore putinterval
  0 1 row1max 1 sub {
   /row ED
   0 1 col2max 1 sub {
    /col ED
    /sum 0 def
    0 1 col1max 1 sub{
    /rowcol ED
    sum
    m1 row col1max mul rowcol add get
    m2 rowcol col2max mul col add get
    mul add 
    /sum ED
    } for
    m3 row col2max mul col add sum put
   } for
  } for
  m3
  end % end of M@tMulDict

} def
%
/SetMQuaternion {

  /MnewTOold 11 array def

  /Qu@ternionDict 30 dict def
  Qu@ternionDict begin

  /normRotVec  xRotVec yRotVec zRotVec 3 array astore VecNorm  def
  normRotVec 0 gt
  {/xRotVecNorm xRotVec normRotVec div def
   /yRotVecNorm yRotVec normRotVec div def
   /zRotVecNorm zRotVec normRotVec div def
   RotAngle}
  {/xRotVecNorm 1 def
   /yRotVecNorm 0 def
   /zRotVecNorm 0 def 
   0} ifelse

  2 div dup
  /q0 exch cos def
      sin dup dup
  /q1 exch xRotVecNorm mul def
  /q2 exch yRotVecNorm mul def
  /q3 exch zRotVecNorm mul def

  /q0q0 q0 q0 mul def
  /q0q1 q0 q1 mul def
  /q0q2 q0 q2 mul def
  /q0q3 q0 q3 mul def

  /q1q1 q1 q1 mul def
  /q1q2 q1 q2 mul def
  /q1q3 q1 q3 mul def

  /q2q2 q2 q2 mul def
  /q2q3 q2 q3 mul def

  /q3q3 q3 q3 mul def

  MnewTOold 0 q0q0 q1q1 add q2q2 sub q3q3 sub put
  MnewTOold 1 q1q2 q0q3 sub 2 mul put
  MnewTOold 2 q1q3 q0q2 add 2 mul put

  MnewTOold 3 q1q2 q0q3 add 2 mul put
  MnewTOold 4 q0q0 q1q1 sub q2q2 add q3q3 sub put
  MnewTOold 5 q2q3 q0q1 sub 2 mul put

  MnewTOold 6 q1q3 q0q2 sub 2 mul put
  MnewTOold 7 q2q3 q0q1 add 2 mul put
  MnewTOold 8 q0q0 q1q1 sub q2q2 sub q3q3 add put

  MnewTOold 9 3 put
  MnewTOold 10 3 put

  end % end of Qu@ternionDict

} def
%
/SetMxyz {
  1.0 0.0 0.0  0.0 1.0 0.0  0.0 0.0 1.0  3 3  11 array astore /MnewTOold ED
  RotSequence cvx exec % Now create a new MnewTOold using xyz, etc.
} def
%
/ConcatMQuaternion {
  MnewTOold % Push onto stack
  SetMQuaternion % Uses [xyz]RotVec and RotAngle to make MnewToOld 
  MnewTOold matmul /MnewTOold ED
} def
%
/ConcatMxyz {
  MnewTOold % Push onto stack
  SetMxyz % Uses RotX, etc. to set MnewTOold 
  MnewTOold matmul /MnewTOold ED
} def
%
/RotatePoint{
  MnewTOold x y z  3 1  5 array astore matmul
  0 3 getinterval aload pop 
  /z ED 
  /y ED 
  /x ED 
} def
%
/makeMoldTOnew {
  /MoldTOnew 11 array def
  MoldTOnew 0 MnewTOold 0 get put
  MoldTOnew 1 MnewTOold 3 get put
  MoldTOnew 2 MnewTOold 6 get put
  MoldTOnew 3 MnewTOold 1 get put
  MoldTOnew 4 MnewTOold 4 get put
  MoldTOnew 5 MnewTOold 7 get put
  MoldTOnew 6 MnewTOold 2 get put
  MoldTOnew 7 MnewTOold 5 get put
  MoldTOnew 8 MnewTOold 8 get put
  MoldTOnew 9               3 put
  MoldTOnew 10              3 put
} def
%
/RotXaxis { 
  eulerRotation 
  {1 0 0}
  {makeMoldTOnew MoldTOnew  1 0 0  3 1  5 array astore matmul
   0 3 getinterval aload pop} ifelse
  /zRotVec ED
  /yRotVec ED
  /xRotVec ED
  /RotAngle RotX def
  ConcatMQuaternion
} def
/RotYaxis { 
  eulerRotation 
  {0 1 0}
  {makeMoldTOnew MoldTOnew  0 1 0  3 1  5 array astore matmul
   0 3 getinterval aload pop} ifelse
  /zRotVec ED
  /yRotVec ED
  /xRotVec ED
  /RotAngle RotY def
  ConcatMQuaternion
} def
/RotZaxis { 
  eulerRotation 
  {0 0 1}
  {makeMoldTOnew MoldTOnew  0 0 1  3 1  5 array astore matmul
   0 3 getinterval aload pop} ifelse
  /zRotVec ED
  /yRotVec ED
  /xRotVec ED
  /RotAngle RotZ def
  ConcatMQuaternion
} def
/xyz { RotXaxis RotYaxis RotZaxis } def
/yxz { RotYaxis RotXaxis RotZaxis } def
/yzx { RotYaxis RotZaxis RotXaxis } def
/xzy { RotXaxis RotZaxis RotYaxis } def
/zxy { RotZaxis RotXaxis RotYaxis } def
/zyx { RotZaxis RotYaxis RotXaxis } def
/quaternion { } def % Null
%
/VecNorm { 0 exch { dup mul add } forall sqrt } def
%
/UnitVec {			% on stack is [a]; returns a vector with [a][a]/|a|=1 
  dup VecNorm /norm ED
  norm 0 lt {/norm 0 def} if
  { norm div } forall 3 array astore } def
%
/AxB {				% on the stack are the two vectors [a][b]
    aload pop /b3 ED /b2 ED /b1 ED
    aload pop /a3 ED /a2 ED /a1 ED
    a2 b3 mul a3 b2 mul sub
    a3 b1 mul a1 b3 mul sub
    a1 b2 mul a2 b1 mul sub
    3 array astore } def
%
/AaddB {			% on the stack are the two vectors [a][b]
    aload pop /b3 ED /b2 ED /b1 ED
    aload pop /a3 ED /a2 ED /a1 ED
    a1 b1 add a2 b2 add a3 b3 add
    3 array astore } def
%
/AmulC {			% on stack is [a] and c; returns [a] mul c
    /factor ED { factor mul } forall 3 array astore } def
%
%
/setColorLight { % expects 7 values on stack C M Y K xL yL zL
% les rayons de lumi�re
  xLight dup mul yLight dup mul zLight dup mul add add sqrt /NormeLight ED
% the color values
  /K ED
  /Yellow ED
  /Magenta ED
  /Cyan ED
} def
%
/facetteSphere {
  newpath
  /Xpoint Rsphere theta cos mul phi cos mul CX add def
  /Ypoint Rsphere theta sin mul phi cos mul CY add def
  /Zpoint Rsphere phi sin mul CZ add def
  Xpoint Ypoint Zpoint tx@3Ddict begin ProjThreeD end moveto
  theta 1 theta increment add {%
    /theta1 ED
    /Xpoint Rsphere theta1 cos mul phi cos mul CX add def
    /Ypoint Rsphere theta1 sin mul phi cos mul CY add def
    /Zpoint Rsphere phi sin mul CZ add def
    Xpoint Ypoint Zpoint tx@3Ddict begin ProjThreeD end  lineto
  } for
  phi 1 phi increment add {
    /phi1 ED
    /Xpoint Rsphere theta increment add cos mul phi1 cos mul CX add def
    /Ypoint Rsphere theta increment add sin mul phi1 cos mul CY add def
    /Zpoint Rsphere phi1 sin mul CZ add def
    Xpoint Ypoint Zpoint tx@3Ddict begin ProjThreeD end lineto
  } for
  theta increment add -1 theta {%
    /theta1 ED
    /Xpoint Rsphere theta1 cos mul phi increment add cos mul CX add def
    /Ypoint Rsphere theta1 sin mul phi increment add cos mul CY add def
    /Zpoint Rsphere phi increment add sin mul CZ add def
    Xpoint Ypoint Zpoint tx@3Ddict begin ProjThreeD end lineto
  } for
  phi increment add -1 phi {
    /phi1 ED
    /Xpoint Rsphere theta cos mul phi1 cos mul CX add def
    /Ypoint Rsphere theta sin mul phi1 cos mul CY add def
    /Zpoint Rsphere phi1 sin mul CZ add def
    Xpoint Ypoint Zpoint tx@3Ddict begin ProjThreeD end lineto
  } for
  closepath 
} def
%
/MaillageSphere { 
% on stack must be x y z Radius increment C M Y K 
  setColorLight
  /increment ED
  /Rsphere ED
  /CZ ED
  /CY ED
  /CX ED
  /StartTheta 0 def
  /condition { PSfacetteSphere 0 ge } def
  -90 increment 90 increment sub {%
    /phi ED
    StartTheta increment 360 StartTheta add increment sub {%
      /theta ED
      % Centre de la facette
      /Xpoint Rsphere theta increment 2 div add cos mul phi increment 2 div add cos mul CX add def
      /Ypoint Rsphere theta increment 2 div add sin mul phi increment 2 div add cos mul CY add def
      /Zpoint Rsphere phi increment 2 div add sin mul CZ add def
      % normale a la facette
      /nXfacette Xpoint CX sub def
      /nYfacette Ypoint CY sub def
      /nZfacette Zpoint CZ sub def
      % test de visibilite
      /PSfacetteSphere 
        vX nXfacette mul
        vY nYfacette mul add
        vZ nZfacette mul add
      def
      condition {
        gsave
        facetteSphere
        /cosV { 1 xLight nXfacette mul
          yLight nYfacette mul
          zLight nZfacette mul
          add add
          NormeLight
          nXfacette dup mul
          nYfacette dup mul
          nZfacette dup mul
          add add sqrt mul div sub } bind def
        Cyan cosV mul Magenta cosV mul Yellow cosV mul K cosV mul setcmykcolor fill 
	grestore
%	0 setgray
        showgrid { facetteSphere stroke } if
      } if 
    } for
    % /StartTheta StartTheta increment 2 div add def
  } for
} def
%
%---------------------- Cylinder ---------------------------
%
/PlanCoupeCylinder { %
  /TableauxPoints [
    0 1 359 { 
      /phi ED 
      [ Radius phi Height ConvCyl2d ] % on décrit le cercle
    } for
  ] def
  newpath
  TableauxPoints 0 get aload pop moveto
  1 1 359 { TableauxPoints exch get aload pop lineto } for
  closepath
} def
%
/facetteCylinder { % 
    newpath
    Radius phi currentHeight ConvCyl2d moveto
    phi 1 phi dAngle add  { % loop variable on stack
      Radius exch currentHeight ConvCyl2d lineto        
    } for
    phi dAngle add -1 phi { %	fill dHeight
      Radius exch currentHeight dHeight add ConvCyl2d lineto 
    } for
    closepath
  } def % facette
%
/MaillageCylinder { % on stack true or false for saving values
    { setColorLight  % expects 4 values on stack C M Y K
      /dHeight ED /dAngle ED /Height ED /Radius ED
      /CZ ED /CY ED /CX ED } if
%     
    0 dHeight Height dHeight sub {
      /currentHeight ED
      0 dAngle 360 dAngle sub {
        /phi ED
% Normal vector of the center
        /nXfacetteCylinder Radius phi dAngle 2 div add cos mul CX add def 
        /nYfacetteCylinder Radius phi dAngle 2 div add sin mul CY add def 
        /nZfacetteCylinder currentHeight dHeight 2 div add CZ add def 
        /NormeN 
          nXfacetteCylinder dup mul
          nYfacetteCylinder dup mul
          nZfacetteCylinder dup mul
          add add sqrt def
        NormeN 0 eq { /NormeN 1e-10 def } if
% test de visibilité
       /PSfacetteCylinder 
    	    vX nXfacetteCylinder mul
            vY nYfacetteCylinder mul add
            vZ nZfacetteCylinder mul add def
       condition {
         facetteCylinder
         /cosV 
	   1 xLight nXfacetteCylinder mul
           yLight nYfacetteCylinder mul
           zLight nZfacetteCylinder mul
           add add
	   NormeLight NormeN mul div sub def
         Cyan Magenta Yellow K
         cosV mul 4 1 roll cosV mul 4 1 roll 
	 cosV dup mul mul 4 1 roll cosV dup mul mul 4 1 roll
         setcmykcolor fill
          showgrid { 
            0 setgray
            facetteCylinder % drawing the segments
            stroke } if
       } if
     } for
    } for
} def
%
%------------------------ Cylinder type II -----------------------
%
/MoveTo { Conv3D2D moveto } def
/LineTo { Conv3D2D lineto } def

/IIIDEllipse { % x y z rA rB startAngle endAngle Wedge
  /dAngle 1 def
  /isWedge ED
  /endAngle ED
  /startAngle ED
  /radiusB ED
  /radiusA ED
  startAngle cos radiusA mul startAngle sin radiusB mul 0 
  isWedge { 0 0 moveto LineTo }{ MoveTo } ifelse
  /Angle startAngle def
  startAngle dAngle endAngle {
    /Angle ED
    Angle cos radiusA mul Angle sin radiusB mul 0 LineTo  
  } for
  isWedge { 0 0 lineto } if
} def

/IIIDCircle { % x y z r startAngle endAngle Wedge
  7 3 roll % startAngle endAngle Wedge x y z r
  dup      % startAngle endAngle Wedge x y z r r
  8 -3 roll
  IIIDEllipse 
} def

/IIIDWedge { % x y z r startAngle endAngle
  true IIIDCircle
} def

/IIIDCylinder {% x y z r h start end wedge
  /isWedge ED
  /increment ED
  /endAngle ED
  /startAngle ED
  /height ED
  /radius ED
  startAngle increment endAngle {
    /Angle ED
    radius Angle 0 ConvCylToCartesian MoveTo  
    radius Angle height ConvCylToCartesian LineTo  
  } for
  stroke
} def
%
%---------------------- Box ---------------------------
%
/PlanCoupeBox { % x y z
  /TableauxPoints [
      [ CX CY CZ Height add ConvBox2d ] % top or bottom
      [ CX CY Depth add CZ Height add ConvBox2d ]
      [ CX Width add CY Depth add CZ Height add ConvBox2d ] 
      [ CX Width add CY CZ Height add ConvBox2d ] 
      [ CX CY CZ Height add ConvBox2d ] % bottom
    ] def
    newpath
    TableauxPoints 0 get aload pop moveto
    0 1 3 {
      TableauxPoints exch get aload pop
      lineto } for
    closepath
} def
%
/facetteBox { % 
    newpath
    dup
    1 eq { % back
      CX CY CZ ConvBox2d moveto
      CX CY CZ Height add ConvBox2d lineto
      CX Width add CY CZ Height add ConvBox2d lineto
      CX Width add CY CZ ConvBox2d lineto
      CX CY CZ ConvBox2d lineto
    } if
    dup
    2 eq { % right
      CX CY CZ ConvBox2d moveto
      CX CY CZ Height add ConvBox2d lineto
      CX CY Depth add CZ Height add ConvBox2d lineto
      CX CY Depth add CZ ConvBox2d lineto
      CX CY CZ ConvBox2d lineto
    } if
    dup
    3 eq { % left
      CX Width add CY CZ ConvBox2d moveto
      CX Width add CY Depth add CZ ConvBox2d lineto
      CX Width add CY Depth add CZ Height add ConvBox2d lineto
      CX Width add CY CZ Height add ConvBox2d lineto
      CX Width add CY CZ ConvBox2d lineto
    } if
    4 eq { % front
      CX CY Depth add CZ ConvBox2d moveto
      CX CY Depth add CZ Height add ConvBox2d lineto
      CX Width add CY Depth add CZ Height add ConvBox2d lineto
      CX Width add CY Depth add CZ ConvBox2d lineto
      CX CY Depth add CZ ConvBox2d lineto
    } if
    closepath
  } def % facette
%
/TestPlane { % on stack x y z of the plane center and # of plane
  /nZfacetteBox ED /nYfacetteBox ED /nXfacetteBox ED
  /Plane ED
  /NormeN 
    nXfacetteBox dup mul
    nYfacetteBox dup mul
    nZfacetteBox dup mul
    add add sqrt def
  NormeN 0 eq { /NormeN 1e-10 def } if
% test de visibilite
  /PSfacetteBox 
    vX nXfacetteBox mul
    vY nYfacetteBox mul add
    vZ nZfacetteBox mul add def
  condition {
    Plane facetteBox
         /cosV 
	   1 xLight nXfacetteBox mul
           yLight nYfacetteBox mul
           zLight nZfacetteBox mul
           add add
	   NormeLight NormeN mul div sub def
         Cyan Magenta Yellow K
         cosV mul 4 1 roll cosV mul 4 1 roll 
	 cosV dup mul mul 4 1 roll cosV dup mul mul 4 1 roll
         setcmykcolor fill
         0 setgray
         Plane facetteBox % drawing the segments
         stroke
       } if
} def
%
/MaillageBox { % on stack true or false for saving values
    { setColorLight  % expects 4 values on stack C M Y K 
      /Depth ED /Height ED /Width ED
      /CZ ED /CY ED /CX ED } if
%
% Normal vector of the box center
  /PlaneSet [
    [ Width 2 div CX add 
      CY 
      Height 2 div CZ add ] % normal back
    [ CX 
      Depth 2 div CY add 
      Height 2 div CZ add ] % normal right
    [ Width CX add 
      Depth 2 div CY add 
      Height 2 div CZ add ] % normal left
    [ Width 2 div CX add 
      Depth CY add 
      Height 2 div CZ add ] % normal front
  ] def
  PlaneSequence length 0 eq { % user defined?
    Alpha abs cvi 360 mod /iAlpha ED
    iAlpha 90 lt { [ 1 2 3 4 ]  
      }{ iAlpha 180 lt { [ 2 4 1 3 ]  
        }{ iAlpha 270 lt { [ 3 4 1 2 ] }{ [ 3 1 4 2] } ifelse } ifelse } ifelse 
  }{ PlaneSequence } ifelse 
  { dup 1 sub PlaneSet exch get aload pop TestPlane } forall
} def
%
%--------------------------- Paraboloid -----------------------------
/PlanCoupeParaboloid {
    /Z height store
    /V {Z sqrt} bind def
    /TableauxPoints [
      0 1 359 { 
        /U ED [ U U Z V calculate2DPoint ] % on decrit le cercle
      } for
    ] def
    newpath
    TableauxPoints 0 get aload pop moveto
    0 1 359 {
      /compteur ED
      TableauxPoints compteur get aload pop
      lineto } for
    closepath
} def
%
/facetteParaboloid{
    newpath
    U U Z V calculate2DPoint moveto
    U 1 U increment add  {%
      /U1 ED
      U1 U1 Z V calculate2DPoint lineto
    } for
    Z pas10 Z pas add pas10 add{
      /Z1 ED
      /V {Z1 sqrt} bind def
      U1 U1 Z1 V calculate2DPoint lineto
    } for
    U increment add -1 U {%
      /U2 ED
      U2 U2 Z pas add V calculate2DPoint lineto
    } for
    Z pas add pas10 sub pas10 neg Z pas10 sub {
      /Z2 ED
      /V Z2 abs sqrt def
      U U Z2 V calculate2DPoint lineto
    } for
    closepath
} def % facette
%
/MaillageParaboloid {
  % on stack true or false for saving values
    { setColorLight  % expects 7 values on stack C M Y K xL yL zL 
%      /CZ ED /CY ED /CX ED 
    } if    
    0 pas height pas sub {%
      /Z ED
      /V Z sqrt def
      0 increment 360 increment sub {%
        /U ED
% Centre de la facette
        /Ucentre U increment 2 div add def
        /Vcentre Z pas 2 div add sqrt def
% normale à la facette
        /nXfacetteParaboloid 2 Vcentre dup mul mul Ucentre cos mul radius mul def
        /nYfacetteParaboloid 2 Vcentre dup mul mul Ucentre sin mul radius mul def
        /nZfacetteParaboloid Vcentre neg radius dup mul mul def
        /NormeN {
          nXfacetteParaboloid dup mul
          nYfacetteParaboloid dup mul
          nZfacetteParaboloid dup mul
          add add sqrt} bind def
        NormeN 0 eq {/NormeN 1e-10 def} if
% test de visibilit�
       /PSfacetteParaboloid vX nXfacetteParaboloid mul
                  vY nYfacetteParaboloid mul add
                  vZ nZfacetteParaboloid mul add def
       condition {
         facetteParaboloid
         /cosV 1 xLight nXfacetteParaboloid mul
           yLight nYfacetteParaboloid mul
           zLight nZfacetteParaboloid mul
           add add
           NormeLight
           NormeN mul div sub def
         Cyan Magenta Yellow K  
         cosV mul 4 1 roll cosV mul 4 1 roll cosV dup mul mul 4 1 roll cosV dup mul mul 4 1 roll
         setcmykcolor fill
         showgrid {
           0 setgray
           facetteParaboloid
           stroke } if
       } if
     } for
    } for
} def
%
% ------------------------------------ math stuff ----------------------------------
%
% Matrix A in arrays of rows A[[row1][row2]...]
% with [row1]=[a11 a12 ... b1]
% returns on stack solution vector X=[x1 x2 ... xn]
/SolveLinEqSystem { 				% on stack matrix M=[A,b] (A*x=b)
  10 dict begin					% hold all ocal
    /A exch def
    /Rows A length def         			% Rows = number of rows
    /Cols A 0 get length def   			% Cols = number of columns
    /Index [ 0 1 Rows 1 sub { } for ] def	% Index = [0 1 2 ... Rows-1]
    /col 0 def
    /row  0 def
    /PR Rows array def 				% PR[c] = pivot row for row row
  { 						% starts the loop, find pivot entry in row r
    col Cols ge row  Rows ge or { exit } if	% col < Cols and row < Rows else exit
    /pRow row def  				% pRow = pivot row		
    /max A row  get col get abs def		% get A[row[col]], first A[0,0] 
    row 1 add 1 Rows 1 sub { 			% starts for loop 1 1 Rows-1
      /j exch def				% index counter
      /x A j get col get abs def		% get A[j[r]]
      x max gt {				% x>max, then save position
        /pRow j def
        /max x def
      } if
    } for					% now we have the row with biggest A[0,1]
						% with pRow = the pivot row
    max 0 gt {					% swap entries pRow and row  in i 
      /tmp Index row  get def
      Index row  Index pRow get put
      Index pRow tmp put			% and columns pRow and row  in A
      /tmp A row get def
      A row  A pRow get put
      A pRow tmp put   				% pivot
      /row0  A row  get def 			% the pivoting row
      /p0 row0  col get def 			% the pivot value
      row 1 add 1 Rows 1 sub { 			% start for loop
        /j exch def
        /c1 A j get def
        /p c1 col get p0 div def
        c1 col p put				% subtract (p1/p0)*row[i] from row[j]
        col 1 add 1 Cols 1 sub {		% start for loop
          /i exch def
          c1 dup i exch 			% c1 i c1
          i get row0 i get p mul sub put
        } for
      } for
      PR row col put
      /col col 1 add def
      /row row 1 add def
    }{						% all zero entries
      /row row 1 add def			% continue loop with same row
    } ifelse
  } loop
  /X A def					% solution vector
  A Rows 1 sub get dup
  Cols 1 sub get exch
  Cols 2 sub get div
  X Rows 1 sub 3 -1 roll put  			% X[n]
  Rows 2 sub -1 0 {				% for loop to calculate X[i]
    /xi exch def				% current index
    A xi get 					% i-th row
    /Axi exch def
    /sum 0 def
    Cols 2 sub -1 xi 1 add { 
      /n exch def
      /sum sum Axi n get X n get mul add def 
    } for
    Axi Cols 1 sub get 				% b=Axi[Cols-1]
    sum sub 					% b-sum
    Axi xi get div				% b-sum / Axi[xi]
    X xi 3 -1 roll put  			% X[xi]
  } for
  X
  end 
} def
%
% u -> e_u with |e_u|=1 
/vector-unit { 1 dict begin
  dup vector-length 1 exch div 
  vector-scale
  end 
} def
%
% u v -> u+v
/vector-add { 1 dict begin
  /v exch def
  [ exch
  0 	     	% u i
  exch { 	% i u[i]
    v 		% i u[i] v
    2 index get add 	% i u[i]+v[i]
    exch 1 add	% i
  } forall
  pop
  ]
  end 
} def
%
% u v -> u-v
/vector-sub { 1 dict begin
  /v exch def
  [ exch
  0 	     	% u i
  exch {	% i u[i]
    v 		% i u[i] v
    2 index get sub 	% i u[i]+v[i]
    exch 1 add	% i
  } forall
  pop
  ]
end } def
%
% [v] c -> [c.v]
/vector-scale { 1 dict begin
  /c exch def
  [ exch
  { 		% s i u[i]
    c mul	% s i u[i] v 
  } forall
  ]
  end } def
%
%
% [u] [v] -> [u x v]
/vector-prod { %% x1 y1 z1 x2 y2 z2
6 dict begin
  aload pop 
  /zp exch def /yp exch def /xp exch def
  aload pop 
  /z exch def /y exch def /x exch def
  [ y zp mul z yp mul sub
   z xp mul x zp mul sub
   x yp mul y xp mul sub ]
end
} def
%
% [u] [v] -> u.v
/vector-mul { %% x1 y1 z1 x2 y2 z2
6 dict begin
  aload pop 
  /zp exch def /yp exch def /xp exch def
  aload pop 
  /z exch def /y exch def /x exch def
  x xp mul y yp mul add z zp mul add
end
} def
%
% [x y z ... ] -> r
% watch out for overflow
/vector-length { 1 dict begin
dup
% find maximum entry
/max 0 def
{ % max 
  abs dup max gt {
    % if abs gt max
    /max exch def
  } {
    pop
  } ifelse
} forall
max 0 ne {
  0 exch 
  {  % 0 v[i]
    max div dup mul add
  } forall
  sqrt
  max mul
} {
  pop 0
} ifelse
end } def
%
end % tx@3DPlotDict
%