home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Transactor
/
Transactor_21_1988_Transactor_Publishing.d64
/
hologram
(
.txt
)
< prev
next >
Wrap
Commodore BASIC
|
2023-02-26
|
10KB
|
375 lines
10 rem this program makes binary holograms
15 rem by: patrick hawley, port elgin, on
20 dim x(32),y(32),s(13),u(13),ar(32,32),ai(32,32),ch$(11,15),e(16)
30 ti$="000000":open4,4
40 ir=1: ii=1
50 read n,it
60 print"[147]"
70 input"do you want a hologram plot (y/n)";hp$
80 input"random or constant phase (r/c)";ip$
90 pq=13
100 for d=1 to n
110 for e=1 to n
120 read ar(d,e)
130 ai(d,e)=ar(d,e)*ii
140 ar(d,e)=ar(d,e)*ir
150 ifip$="c"then200
160 ar(d,e)=(-1+2*rnd(1))*sqr(ir^2+ii^2)*ar(d,e)
170 pm=int(-1+3*rnd(1))
180 if pm=0 then 170
190 ai(d,e)=sqr(ir^2+ii^2-ar(d,e)^2)*ai(d,e)*pm
200 next e,d
220 data 32,1
230 data 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
240 data 0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0
250 data 0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0
260 data 0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0
270 data 0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0
280 data 0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0
290 data 0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0
300 data 0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0
310 data 0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0
320 data 0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0
330 data 0,0,1,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,1,1,0,0
340 data 0,0,1,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,1,1,0,0
350 data 0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0
360 data 0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0
370 data 0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0
380 data 0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0
390 data 0,1,1,0,0,0,1,1,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,1,0,0,0,1,1,0
400 data 0,1,1,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,1,1,0
410 data 0,1,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,1,0
420 data 0,0,1,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,1,0,0
430 data 0,0,1,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,1,0,0
440 data 0,0,1,1,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,1,1,0,0
450 data 0,0,0,1,1,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,1,1,0,0,0
460 data 0,0,0,1,1,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,1,1,0,0,0
470 data 0,0,0,0,1,1,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,1,1,0,0,0,0
480 data 0,0,0,0,0,1,1,0,0,0,0,0,0,1,1,1,1,1,1,0,0,0,0,0,0,1,1,0,0,0,0,0
490 data 0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0
500 data 0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0
510 data 0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0
520 data 0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0
530 data 0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0
540 data 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
550 :
560 gosub 660:rem fast fourier transform
570 gosub 2130:rem amplitude and angle
580 gosub 2370:rem normalize/truncate
590 gosub 2900:rem hologram plot
600 gosub 4340:rem find components for inverse fft from plot data
610 it=-1
620 gosub 660:rem inverse fft
630 gosub 4430:rem gray-scale plot
640 print"finished. time is ";ti$
650 end
660 rem 2d fft
670 for rw=1 to n
680 print" working on fft;row:";rw
690 for re=1 to n
700 x(re)=ar(rw,re):y(re)=ai(rw,re)
710 next re
720 gosub 890
730 for re=1 to n
740 ar(rw,re)=x(re):ai(rw,re)=y(re)
750 next re,rw
770 for cl=1 to n
780 print" working on fft;column:";cl
790 for ce=1 to n
800 x(ce)=ar(ce,cl):y(ce)=ai(ce,cl)
810 next ce
820 gosub 890
830 for ce=1 to n
840 ar(ce,cl)=x(ce):ai(ce,cl)=y(ce)
850 next ce,cl
870 return
880 :
890 rem complex fast fourier transform
900 :
910 if it>0 then 940
920 for i=1 to n
930 y(i)=-y(i): next i
940 lg=log(n)/log(2)
950 if lg<=1 then 1390
960 for k=2 to lg step 2
970 m=2^(lg-k)
980 m4=4*m
990 for j=1 to m
1000 ar=(j-1)/m4*2*(NULL)
1010 c1=cos(ar)
1020 s1=sin(ar)
1030 c2=c1*c1-s1*s1
1040 s2=c1*s1+c1*s1
1050 c3=c2*c1-s2*s1
1060 s3=c2*s1+s2*c1
1070 for i=m4 to n step m4
1080 j0=i+j-m4
1090 j1=j0+m
1100 j2=j1+m
1110 j3=j2+m
1120 r0=x(j0)+x(j2)
1130 r1=x(j0)-x(j2)
1140 i0=y(j0)+y(j2)
1150 i1=y(j0)-y(j2)
1160 r2=x(j1)+x(j3)
1170 r3=x(j1)-x(j3)
1180 i2=y(j1)+y(j3)
1190 i3=y(j1)-y(j3)
1200 x(j0)=r0+r2
1210 y(j0)=i0+i2
1220 if ar=0 then 1300
1230 x(j2)=(r1+i3)*c1+(i1-r3)*s1
1240 y(j2)=(i1-r3)*c1-(r1+i3)*s1
1250 x(j1)=(r0-r2)*c2+(i0-i2)*s2
1260 y(j1)=(i0-i2)*c2-(r0-r2)*s2
1270 x(j3)=(r1-i3)*c3+(i1+r3)*s3
1280 y(j3)=(i1+r3)*c3-(r1-i3)*s3
1290 goto 1360
1300 x(j2)=r1+i3
1310 y(j2)=i1-r3
1320 x(j1)=r0-r2
1330 y(j1)=i0-i2
1340 x(j3)=r1-i3
1350 y(j3)=i1+r3
1360 next i,j,k
1390 if lg=int(lg/2)*2 then 1500
1400 for i=1 to n step 2
1410 r0=x(i)+x(i+1)
1420 r1=x(i)-x(i+1)
1430 i0=y(i)+y(i+1)
1440 i1=y(i)-y(i+1)
1450 x(i)=r0
1460 y(i)=i0
1470 x(i+1)=r1
1480 y(i+1)=i1
1490 next i
1500 s(13)=n/2
1510 u(13)=n
1520 for k=2 to 12
1530 j=14-k
1540 s(j)=1
1550 u(j)=s(j+1)
1560 if s(j+1)>1 then s(j)=int(s(j+1)/2)
1570 next k
1580 al=s(2)
1590 jj=0
1600 a=1: b=a: c=b: d=c: e=d: f=e
1660 for g=f to u(7) step s(7)
1670 for h=g to u(8) step s(8)
1680 for i=h to u(9) step s(9)
1690 for j=i to u(10) step s(10)
1700 for k=j to u(11) step s(11)
1710 for l=k to u(12) step s(12)
1720 for m=l to u(13) step s(13)
1730 jj=jj+1
1740 if jj<=m then 1810
1750 t=x(jj)
1760 x(jj)=x(m)
1770 x(m)=t
1780 t=y(jj)
1790 y(jj)=y(m)
1800 y(m)=t
1810 :
1820 next m,l,k,j,i,h,g
1890 f=f+s(6)
1900 if f<=u(6) then 1660
1910 e=e+s(5)
1920 if e<=u(5) then 1650
1930 d=d+s(4)
1940 if d<=u(4) then 1640
1950 c=c+s(3)
1960 if c<=u(s) then 1630
1970 b=b+s(2)
1980 if b<=u(2) then 1620
1990 a=a+1
2000 if a<=al then 1610
2010 for sl=1 to (n-2)/2
2020 bc=n+1-sl:fc=sl+1
2030 tx=x(fc):ty=y(fc)
2040 x(fc)=x(bc):y(fc)=y(bc)
2050 x(bc)=tx:y(bc)=ty
2060 next sl
2070 if it>0 then return
2080 for i=1 to n
2090 x(i)=x(i)/n
2100 y(i)=-y(i)/n
2110 nexti
2120 return
2130 rem convert into amplitude & angle
2140 print"...finding amplitude & angle"
2150 for d=1 to n
2160 for e=1 to n
2170 ap=sqr(ar(d,e)*ar(d,e)+ai(d,e)*ai(d,e))
2180 if ap=0 then 2330
2190 if ar(d,e)<>0 then 2230
2200 if ai(d,e)<0 then pa=-(NULL)/2
2210 if ai(d,e)>0 then pa=(NULL)/2
2220 goto 2320
2230 if ai(d,e)<>0 then 2270
2240 if ar(d,e)<0 the npa=(NULL)
2250 if ar(d,e)>0 thenpa=0
2260 goto 2320
2270 pa=atn(ai(d,e)/ar(d,e))
2280 if ar(d,e)>0 and ai(d,e)>0 thenpa=pa
2290 if ar(d,e)<0 and ai(d,e)>0 thenpa=(NULL)+pa
2300 if ar(d,e)<0 and ai(d,e)<0 thenpa=pa-(NULL)
2310 if ar(d,e)>0 and ai(d,e)<0 thenpa=pa
2320 ar(d,e)=ap: ai(d,e)=pa
2330 next e,d
2350 return
2360 :
2370 rem amplitude normalization routine
2380 rem find largest amplitude
2390 la=ar(1,1)
2400 for d=1 to n
2410 for e=1 to n
2420 ifar(d,e)>lathenla=ar(d,e)
2430 next e,d
2450 if ip$="r" then 2820
2460 print"the largest amplitude is";la
2470 rem check amplitude distribution
2480 p8=0: p6=0: p4=0: p2=0: p1=0
2490 n8=.8*la:n6=.6*la:n4=.4*la:n2=.2*la:n1=.1*la
2500 print" "
2510 for d=1 to n
2520 for e=1 to n
2530 if ar(d,e)<=n8 then p8=p8+1
2540 if ar(d,e)<=n6 then p6=p6+1
2550 if ar(d,e)<=n4 then p4=p4+1
2560 if ar(d,e)<=n2 then p2=p2+1
2570 if ar(d,e)<=n1 then p1=p1+1
2580 next e,d
2600 p8=int(p8/1024*100):p6=int(p6/1024*100):p4=int(p4/1024*100)
2610 p2=int(p2/1024*100):p1=int(p1/1024*100)
2620 print"80% of the max. amp is";n8
2630 print p8;"% of amp's are equal or smaller"
2640 print" "
2650 print"60% of the max. amp is";n6
2660 printp6;"% of amp's are equal or smaller"
2670 print" "
2680 print"40% of the max. amp is";n4
2690 printp4;"% of amp's are equal or smaller"
2700 print" "
2710 print"20% of the max. amp is";n2
2720 printp2;"% of amp's are equal or smaller"
2730 print" "
2740 print"10% of the max. amp is";n1
2750 printp1;"% of amp's are equal or smaller"
2760 print" "
2770 print"time is ";ti$
2780 print" "
2790 input"what no. do you want to normalize to";nl
2800 print" "
2810 input"which means what % of data truncated";pt
2820 ifip$="r"thennl=la
2830 for d=1 to n
2840 for e=1 to n
2850 ar(d,e)=ar(d,e)/nl: if ar(d,e)>1 then ar(d,e)=1
2860 next e,d
2880 return
2890 :
2900 rem plot routine
2910 rem read in plot characters
2920 if hp$="n" then 3720
2930 print"...defining plot characters"
2940 for i=1 to 8:e(i)=2^i-1:e(i+8)=e(i)*2^(8-i):next
2950 for i=1 to 11: for j=1 to 15: for k=1 to 8
2960 ch=0: if k<(13-i) and k>(8-i) then ch=e(j)
2970 ch$(i,j)=ch$(i,j)+chr$(ch)
2980 next k,j,i
3710 :
3720 open 2,4,5
3730 open 6,4,6:print#6,chr$(21)
3740 rem main plotting loop
3750 tb=0
3760 nf=2: if hp$="n" then nf=1
3770 for d=1 to n
3780 for f=1 to nf
3790 for e=1 to n
3800 rem determine which character
3810 rem find character height
3820 ht=8*ar(d,e)
3830 dh=ht-int(ht)
3840 if dh<.125/2 then ht=int(ht):goto3860
3850 ht=int(ht)+1
3860 if f=nf then ar(d,e)=ht/8*nl
3870 rem find phase
3880 cp=pq*ai(d,e)/(2*(NULL))
3890 qp=int(pq/2)
3900 ifabs(cp)-int(abs(cp))>.5then3950
3910 if cp<0 then cp=int(cp):if f=nf then ai(d,e)=(cp+1)*(NULL)/qp
3920 if cp>=0 then cp=int(cp)+1:if f=nf then ai(d,e)=(cp-1)*(NULL)/qp
3930 if hp$="n"then 4290
3940 goto 3980
3950 if cp<0 then cp=int(cp)-1:if f=nf then ai(d,e)=(cp+1)*(NULL)/qp
3960 if cp>0then cp=int(cp)+2:if f=nfthen ai(d,e)=(cp-1)*(NULL)/qp
3970 if hp$="n" then 4290
3980 if f<>1 then 4230
3990 if ht=0 then print#4," ";chr$(141);:goto4260
4000 if abs(cp)>2then 4130
4010 if cp>-2 then 4050
4020 print#2,ch$(3,ht):print#4,chr$(254);chr$(141);
4030 tb=tb+1:print#2,ch$(11,ht):print#4,tab(tb);chr$(254);chr$(141);
4040 goto 4260
4050 if cp>1then 4090
4060 print#2,ch$(2,ht):print#4,chr$(254);chr$(141);
4070 tb=tb+1:print#2,ch$(10,ht):print#4,tab(tb);chr$(254);chr$(141);
4080 goto 4260
4090 print#2,ch$(1,ht):print#4,chr$(254);chr$(141);
4100 tb=tb+1:print#2,ch$(9,ht):print#4,tab(tb);chr$(254);chr$(141);
4110 tb=tb+1:print#2,ch$(9,ht):print#4,tab(tb);chr$(254);chr$(141);
4120 goto 4260
4130 if cp<0 then 4180
4140 cp=cp+11-cp*2
4150 print#4," ";chr$(141);
4160 tb=tb+1:print#2,ch$(cp,ht):print#4,tab(tb);chr$(254);chr$(141);
4170 goto4260
4180 if cp>0 then 4230
4190 cp=abs(cp)+1
4200 print#2,ch$(cp,ht):print#4,chr$(254);chr$(141);
4210 tb=tb+1:print#4,tab(tb);" ";chr$(141);
4220 goto 4260
4230 if ht=0 then 3990
4240 if ht<8 then ht=ht+8
4250 goto 3990
4260 if e=n then print#4,chr$(141);chr$(13);
4270 tb=2*e:if e=n then tb=0
4280 if e<>n then print#4,tab(tb);
4290 next e,f,d
4320 return
4330 :
4340 rem find real and imaginary components from truncated & quantized data
4350 for d=1 to n: for e=1 to n
4370 ap=ar(d,e)
4380 ar(d,e)=ap*cos(ai(d,e))
4390 ai(d,e)=ap*sin(ai(d,e))
4400 next e,d
4420 return
4430 :
4440 rem gray scale picture routine
4450 print#4," "
4460 for d=1 to n
4470 print#4,tab(23);
4480 for e=1 to n
4490 in=sqr(ar(d,e)*ar(d,e)+ai(d,e)*ai(d,e))
4500 af=1/sqr(ir*ir+ii*ii):in=af*in
4510 if in>1 then in=1
4520 if in<.1 thenprint#4," ";:goto4590
4530 if in<.2 thenprint#4,"'";:goto4590
4540 if in<.4 thenprint#4,":";:goto4590
4550 if in<.6 thenprint#4,"*";:goto4590
4560 if in<.8 thenprint#4,"[166]";:goto4590
4570 if in<.9 thenprint#4,"*[146]";:goto4590
4580 if in<=1 thenprint#4," [146]";
4590 if e=n then print#4,chr$(13);
4600 next e,d
4620 print#4," "
4630 print#4,"input real amplitude:";ir
4640 print#4,"input imaginary amplitude:";ii
4650 print#4,"phase angle:";ip$
4660 print#4,"no. of phase quantization levels:";pq
4670 if ip$="r" then 4700
4680 print#4,"transform normalized with...";nl
4690 print#4,"percent of transformed amp's truncated...";pt
4700 return