PROGRAM ROTACIJA_Z_UPORABO_KVATERNIONOV include 'phigsdef.f' real WindowLimits(4)/-20.0,20.0,-20.0,20.0/ real ViewportLimits(4)/0.0,0.3,0.0,0.3/ real ClipLimits(4) /0.0,0.3,0.0,0.3/ real ViewMappingMatrix(3,3) integer wkid, ErrorReturn, i,j,k,l real Qr(4,4),P(8,4),Pr(8,4),px(2),py(2),pz(2) real e1,e2,e3,a,b,h,vs,s1,s2,s3 real*8 u1,u2,u3,sp,ck,cs,w,theta,x,y,z,th,ce,c,theta1 integer R(12,2) c---- Qr rotacijska matrika c---- P tocke telesa c---- u(3) enotski vektor okoli katerega rotira c---- theta kot zasuka c---- x,y,z,w so elementi kvatriniona c---- Pr koordinate tock po rotaciji c---- stevilo vmesnih stopenj do koncnega polozaja c---- ROTACIJA S KVATERNIONI c---- Podaj vrednosti enotskega smernega vektorja, kota zasuka c in koordinato tocke za rotacijo. c WRITE(*,*)'1.Podaj vrednosti enotskega smernega vek. u=x,y,z!' c READ(*,*)u1,u2,u3 WRITE(*,*)'1.Podaj vred. ZACETNEGA enot. smernega vek. S=x,y,z!' READ(*,*)s1,s2,s3 WRITE(*,*)'2.Podaj vred. KONCNEGA enot. smernega vek. E=x,y,z!' READ(*,*)e1,e2,e3 c WRITE(*,*)'3.Podaj vrednosti kota zasuka theta v stopinjah!' c READ(*,*)theta1 WRITE(*,*)'4.Velikost spodnje stranice a 4-str. 3D lika! Max.100' READ(*,*)a WRITE(*,*)'5.Velikost zgornje stranice b 4-str. 3D lika! Max.100' READ(*,*)b WRITE(*,*)'6. Visina 4-str. 3D lika! Max.100' READ(*,*)h WRITE(*,*)'7. Stevilo vmesnih stopenj do koncnega polozaja!' READ(*,*)vs c---- Koordinate tocke za rotacijo P[xp yp zp 0] P(1,1)=+a/2. P(1,2)=-a/2. P(1,3)=-h/2. P(1,4)=+a/2. P(2,1)=+a/2. P(2,2)=+a/2. P(2,3)=-h/2. P(2,4)=0. P(3,1)=-a/2. P(3,2)=a/2. P(3,3)=-h/2. P(3,4)=0. P(4,1)=-a/2. P(4,2)=-a/2. P(4,3)=-h/2. P(4,4)=0. P(5,1)=b/2. P(5,2)=-b/2. P(5,3)=h*2./4. P(5,4)=0. P(6,1)=b/2. P(6,2)=b/2. P(6,3)=h*2./4. P(6,4)=0. P(7,1)=-b/2. P(7,2)=b/2. P(7,3)=h*2./4. P(7,4)=0. P(8,1)=-b/2. P(8,2)=-b/2. P(8,3)=h*2./4. P(8,4)=0. c---- Robovi popisani s tockami R(1,1)=1 R(1,2)=2 R(2,1)=2 R(2,2)=3 R(3,1)=3 R(3,2)=4 R(4,1)=4 R(4,2)=1 R(5,1)=5 R(5,2)=6 R(6,1)=6 R(6,2)=7 R(7,1)=7 R(7,2)=8 R(8,1)=8 R(8,2)=5 R(9,1)=1 R(9,2)=5 R(10,1)=2 R(10,2)=6 R(11,1)=3 R(11,2)=7 R(12,1)=4 R(12,2)=8 c Odpre PHIGS (PhigsOPenPHigs) call popph(1,0) c Odpre delovno postajo (Phigs OPen Workstation) wkid=1 call popwk(wkid," ",WK211280) call pevmm(WindowLimits, ViewportLimits, ErrorReturn, * ViewMappingMatrix) c Nastavitev uporabniskega koordinatnega sistema c ( PhigsSetWorKstationWindow ) c ( wkid, xmin, xmax, ymin, ymax ) c call pswkw(wkid,0.0,1.0,0.0,1.0) c Nastavitev zaslonskih koordinat c ( PhigsSetWorKstationViewport ) c call pswkv(wkid,0.0,0.2,0.0,0.2) c Dolocitev barve crte c ( PhigsSetPolyLineColorIndex ) c Indeksi: c 1 = bela, 2 = rdeca, 3 = rumena, c 4 = zelena, 5 = svetlo modra, 6 = modra, c 7 = violicna, 8 = crna (ozadje) c================ ALGORITEM ZA IZRIS ========================= c==== Risanje lika pred rotacijo -------------------------------- call psplci(5) call pslwsc(1) call psln(pldot) call pswkw(wkid,0.0,1.0,0.0,1.0) call pschh(0.005) c Naris call pswkv(wkid,0.23,0.33,0.2,0.3) do j=1,12 px(1)=P(R(j,1),1) px(2)=P(R(j,2),1) pz(1)=P(R(j,1),3) pz(2)=P(R(j,2),3) call ppl(2,px,pz) enddo c pause c Tloris call pswkv(wkid,0.23,0.33,0.07,0.17) do j=1,12 px(1)=P(R(j,1),1) px(2)=P(R(j,2),1) py(1)=P(R(j,1),2) py(2)=P(R(j,2),2) call ppl(2,px,py) enddo c pause c Stranski ris call pswkv(wkid,0.1,0.2,0.2,0.3) do j=1,12 py(1)=P(R(j,1),2) py(2)=P(R(j,2),2) pz(1)=P(R(j,1),3) pz(2)=P(R(j,2),3) call ppl(2,py,pz) enddo pause c==== Risanje lika po rotaciji ---------------------------- cs=1/SQRT(s1**2+s2**2+s3**2) s1=s1*cs s2=s2*cs s3=s3*cs ce=1/SQRT(e1**2+e2**2+e3**2) e1=e1*ce e2=e2*ce e3=e3*ce sp=(s1*e1)+(s2*e2)+(s3*e3) c---- Pri podanem zacetnem in koncnem vektorju izracunam rotacijski vektor U u1=(s2*e3)-(s3*e2) u2=(s3*e1)-(s1*e3) u3=(s1*e2)-(s2*e1) if ((s1).EQ.(-e1))then u1=0. u2=0. u3=1. elseif ((s2).EQ.(-e2))then u1=0. u2=0. u3=1. elseif((s3).EQ.(-e3))then u1=1. u2=0. u3=0. elseif ((s1).EQ.(e1))then u1=1. u2=0. u3=0. elseif ((s2).EQ.(e2))then u1=0. u2=1. u3=0. elseif((s3).EQ.(e3))then u1=0. u2=0. u3=1. endif c if (abs(s1).EQ.abs(e1))then c u1=0.0 c u2=1. c u3=0.0 c elseif (abs(s2).EQ.abs(e2))then c u1=1.0 c u2=0.0 c u3=0.0 c elseif (abs(s3).EQ.abs(e3))then c u1=0.0 c u2=0.0 c u3=1. c endif c=1/SQRT(u1**2+u2**2+u3**2) u1=u1*c u2=u2*c u3=u3*c c if (abs(s2).EQ.abs(e2).AND.abs(s3).EQ.abs(e3))then c u1=1. c u2=0.0 c u3=0.0 c elseif (abs(s1).EQ.abs(e1).AND.abs(s3).EQ.abs(e3))then c u1=0.0 c u2=1. c u3=0.0 c elseif (abs(s1).EQ.abs(e1).AND.abs(s2).EQ.abs(e2))then c u1=0.0 c u2=0.0 c u3=1. c endif theta1=acos(sp) ck=theta1/3.14159265358979323846264338328*180. write(*,*)'Celoten kot je',ck c---- Izracunavanje vrednosti elementov q=[x y z w] c Enotski smerni vektor mora imeti dolzino 1 zato sem dodal se faktor c. c th=theta1*3.14159265358979323846264338328/180. th=theta1 call psln(plsoli) do l=1,vs c BRISANJE PREJSNJEGA POLOZAJA call psplci(8) c Naris call pswkv(wkid,0.23,0.33,0.2,0.3) do j=1,12 px(1)=Pr(R(j,1),1) px(2)=Pr(R(j,2),1) pz(1)=Pr(R(j,1),3) pz(2)=Pr(R(j,2),3) call ppl(2,px,pz) enddo c Tloris call pswkv(wkid,0.23,0.33,0.07,0.17) do j=1,12 px(1)=Pr(R(j,1),1) px(2)=Pr(R(j,2),1) py(1)=Pr(R(j,1),2) py(2)=Pr(R(j,2),2) call ppl(2,px,py) enddo c Stranski ris call pswkv(wkid,0.1,0.2,0.2,0.3) do j=1,12 py(1)=Pr(R(j,1),2) py(2)=Pr(R(j,2),2) pz(1)=Pr(R(j,1),3) pz(2)=Pr(R(j,2),3) call ppl(2,py,pz) enddo theta=th/vs*l c=1/SQRT(u1**2+u2**2+u3**2) x=u1*sin(theta/2.)*c y=u2*sin(theta/2.)*c z=u3*sin(theta/2.)*c w=cos(theta/2.) c---- Rotacijska matrika Qr c---- Qr(n,m) n-vrsta, m-kolona Qr(1,1)=w**2.+x**2.-y**2.-z**2. Qr(1,2)=2.*x*y+2.*w*z Qr(1,3)=2.*x*z-2.*w*y Qr(1,4)=0. Qr(2,1)=2.*x*y-2.*w*z Qr(2,2)=w**2.-x**2.+y**2.-z**2. Qr(2,3)=2.*y*z+2.*w*x Qr(2,4)=0. Qr(3,1)=2.*x*z+2.*w*y Qr(3,2)=2.*y*z-2.*w*x Qr(3,3)=w**2.-x**2.-y**2.+z**2. Qr(3,4)=0. Qr(4,1)=0. Qr(4,2)=0. Qr(4,3)=0. Qr(4,4)=w**2.+x**2.+y**2.+z**2. call pslwsc(2.) call psln(plsoli) do i=1,8 do k=1,4 Pr(i,k)=Qr(1,k)*P(i,1)+Qr(2,k)*P(i,2)+Qr(3,k)*P(i,3)+Qr(4,k)*0. enddo enddo c Naris call pswkv(wkid,0.23,0.33,0.2,0.3) call psplci(2) do j=1,12 px(1)=Pr(R(j,1),1) px(2)=Pr(R(j,2),1) pz(1)=Pr(R(j,1),3) pz(2)=Pr(R(j,2),3) call ppl(2,px,pz) enddo call ptx(0.2,0.3,'Naris(ravnina X-Z)') c Tloris call pswkv(wkid,0.23,0.33,0.07,0.17) call psplci(4) do j=1,12 px(1)=Pr(R(j,1),1) px(2)=Pr(R(j,2),1) py(1)=Pr(R(j,1),2) py(2)=Pr(R(j,2),2) call ppl(2,px,py) enddo call ptx(0.1,0.2,'Tloris (ravnina X-Y)') c Stranski ris call pswkv(wkid,0.1,0.2,0.2,0.3) call psplci(3) do j=1,12 py(1)=Pr(R(j,1),2) py(2)=Pr(R(j,2),2) pz(1)=Pr(R(j,1),3) pz(2)=Pr(R(j,2),3) call ppl(2,py,pz) enddo call ptx(0.1,0.3,'Stranski ris (ravnina Y-Z)') c pause enddo pause c Zapre delovno postajo c ( PhigsCLoseWorKstation ) call pclwk(wkid) c Zapre PHIGS (zakljuci z graficnim nacinom dela) c ( PhigsCLosePHigs ) call pclph() STOP END