home *** CD-ROM | disk | FTP | other *** search
/ Transactor / Transactor_21_1988_Transactor_Publishing.d64 / hologram (.txt) < prev    next >
Encoding:
Commodore BASIC  |  1988-01-01  |  9.6 KB  |  375 lines
  1. 10 rem this program makes binary holograms
  2. 15 rem by:  patrick hawley, port elgin, on
  3. 20 dim x(32),y(32),s(13),u(13),ar(32,32),ai(32,32),ch$(11,15),e(16)
  4. 30 ti$="000000":open4,4
  5. 40 ir=1: ii=1
  6. 50 read n,it
  7. 60 print"[147]"
  8. 70 input"do you want a hologram plot (y/n)";hp$
  9. 80 input"random or constant phase (r/c)";ip$
  10. 90 pq=13
  11. 100 for d=1 to n
  12. 110 for e=1 to n
  13. 120 read ar(d,e)
  14. 130 ai(d,e)=ar(d,e)*ii
  15. 140 ar(d,e)=ar(d,e)*ir
  16. 150 ifip$="c"then200
  17. 160 ar(d,e)=(-1+2*rnd(1))*sqr(ir^2+ii^2)*ar(d,e)
  18. 170 pm=int(-1+3*rnd(1))
  19. 180 if pm=0 then 170
  20. 190 ai(d,e)=sqr(ir^2+ii^2-ar(d,e)^2)*ai(d,e)*pm
  21. 200 next e,d
  22. 220 data 32,1
  23. 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
  24. 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
  25. 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
  26. 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
  27. 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
  28. 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
  29. 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
  30. 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
  31. 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
  32. 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
  33. 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
  34. 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
  35. 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
  36. 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
  37. 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
  38. 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
  39. 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
  40. 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
  41. 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
  42. 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
  43. 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
  44. 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
  45. 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
  46. 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
  47. 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
  48. 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
  49. 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
  50. 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
  51. 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
  52. 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
  53. 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
  54. 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
  55. 550 :
  56. 560 gosub 660:rem fast fourier transform
  57. 570 gosub 2130:rem amplitude and angle
  58. 580 gosub 2370:rem normalize/truncate
  59. 590 gosub 2900:rem hologram plot
  60. 600 gosub 4340:rem find components for inverse fft from plot data
  61. 610 it=-1
  62. 620 gosub 660:rem inverse fft
  63. 630 gosub 4430:rem gray-scale plot
  64. 640 print"finished. time is  ";ti$
  65. 650 end
  66. 660 rem 2d fft
  67. 670 for rw=1 to n
  68. 680 print"  working on fft;row:";rw
  69. 690 for re=1 to n
  70. 700 x(re)=ar(rw,re):y(re)=ai(rw,re)
  71. 710 next re
  72. 720 gosub 890
  73. 730 for re=1 to n
  74. 740 ar(rw,re)=x(re):ai(rw,re)=y(re)
  75. 750 next re,rw
  76. 770 for cl=1 to n
  77. 780 print"  working on fft;column:";cl
  78. 790 for ce=1 to n
  79. 800 x(ce)=ar(ce,cl):y(ce)=ai(ce,cl)
  80. 810 next ce
  81. 820 gosub 890
  82. 830 for ce=1 to n
  83. 840 ar(ce,cl)=x(ce):ai(ce,cl)=y(ce)
  84. 850 next ce,cl
  85. 870 return
  86. 880 :
  87. 890 rem  complex fast fourier transform
  88. 900 :
  89. 910 if it>0 then 940
  90. 920 for i=1 to n
  91. 930 y(i)=-y(i): next i
  92. 940 lg=log(n)/log(2)
  93. 950 if lg<=1 then 1390
  94. 960 for k=2 to lg step 2
  95. 970 m=2^(lg-k)
  96. 980 m4=4*m
  97. 990 for j=1 to m
  98. 1000 ar=(j-1)/m4*2*(NULL)
  99. 1010 c1=cos(ar)
  100. 1020 s1=sin(ar)
  101. 1030 c2=c1*c1-s1*s1
  102. 1040 s2=c1*s1+c1*s1
  103. 1050 c3=c2*c1-s2*s1
  104. 1060 s3=c2*s1+s2*c1
  105. 1070 for i=m4 to n step m4
  106. 1080 j0=i+j-m4
  107. 1090 j1=j0+m
  108. 1100 j2=j1+m
  109. 1110 j3=j2+m
  110. 1120 r0=x(j0)+x(j2)
  111. 1130 r1=x(j0)-x(j2)
  112. 1140 i0=y(j0)+y(j2)
  113. 1150 i1=y(j0)-y(j2)
  114. 1160 r2=x(j1)+x(j3)
  115. 1170 r3=x(j1)-x(j3)
  116. 1180 i2=y(j1)+y(j3)
  117. 1190 i3=y(j1)-y(j3)
  118. 1200 x(j0)=r0+r2
  119. 1210 y(j0)=i0+i2
  120. 1220 if ar=0 then 1300
  121. 1230 x(j2)=(r1+i3)*c1+(i1-r3)*s1
  122. 1240 y(j2)=(i1-r3)*c1-(r1+i3)*s1
  123. 1250 x(j1)=(r0-r2)*c2+(i0-i2)*s2
  124. 1260 y(j1)=(i0-i2)*c2-(r0-r2)*s2
  125. 1270 x(j3)=(r1-i3)*c3+(i1+r3)*s3
  126. 1280 y(j3)=(i1+r3)*c3-(r1-i3)*s3
  127. 1290 goto 1360
  128. 1300 x(j2)=r1+i3
  129. 1310 y(j2)=i1-r3
  130. 1320 x(j1)=r0-r2
  131. 1330 y(j1)=i0-i2
  132. 1340 x(j3)=r1-i3
  133. 1350 y(j3)=i1+r3
  134. 1360 next i,j,k
  135. 1390 if lg=int(lg/2)*2 then 1500
  136. 1400 for i=1 to n step 2
  137. 1410 r0=x(i)+x(i+1)
  138. 1420 r1=x(i)-x(i+1)
  139. 1430 i0=y(i)+y(i+1)
  140. 1440 i1=y(i)-y(i+1)
  141. 1450 x(i)=r0
  142. 1460 y(i)=i0
  143. 1470 x(i+1)=r1
  144. 1480 y(i+1)=i1
  145. 1490 next i
  146. 1500 s(13)=n/2
  147. 1510 u(13)=n
  148. 1520 for k=2 to 12
  149. 1530 j=14-k
  150. 1540 s(j)=1
  151. 1550 u(j)=s(j+1)
  152. 1560 if s(j+1)>1 then s(j)=int(s(j+1)/2)
  153. 1570 next k
  154. 1580 al=s(2)
  155. 1590 jj=0
  156. 1600 a=1: b=a: c=b: d=c: e=d: f=e
  157. 1660 for g=f to u(7) step s(7)
  158. 1670 for h=g to u(8) step s(8)
  159. 1680 for i=h to u(9) step s(9)
  160. 1690 for j=i to u(10) step s(10)
  161. 1700 for k=j to u(11) step s(11)
  162. 1710 for l=k to u(12) step s(12)
  163. 1720 for m=l to u(13) step s(13)
  164. 1730 jj=jj+1
  165. 1740 if jj<=m then 1810
  166. 1750 t=x(jj)
  167. 1760 x(jj)=x(m)
  168. 1770 x(m)=t
  169. 1780 t=y(jj)
  170. 1790 y(jj)=y(m)
  171. 1800 y(m)=t
  172. 1810 :
  173. 1820 next m,l,k,j,i,h,g
  174. 1890 f=f+s(6)
  175. 1900 if f<=u(6) then 1660
  176. 1910 e=e+s(5)
  177. 1920 if e<=u(5) then 1650
  178. 1930 d=d+s(4)
  179. 1940 if d<=u(4) then 1640
  180. 1950 c=c+s(3)
  181. 1960 if c<=u(s) then 1630
  182. 1970 b=b+s(2)
  183. 1980 if b<=u(2) then 1620
  184. 1990 a=a+1
  185. 2000 if a<=al then 1610
  186. 2010 for sl=1 to (n-2)/2
  187. 2020 bc=n+1-sl:fc=sl+1
  188. 2030 tx=x(fc):ty=y(fc)
  189. 2040 x(fc)=x(bc):y(fc)=y(bc)
  190. 2050 x(bc)=tx:y(bc)=ty
  191. 2060 next sl
  192. 2070 if it>0 then return
  193. 2080 for i=1 to n
  194. 2090 x(i)=x(i)/n
  195. 2100 y(i)=-y(i)/n
  196. 2110 nexti
  197. 2120 return
  198. 2130 rem convert into amplitude & angle
  199. 2140 print"...finding amplitude & angle"
  200. 2150 for d=1 to n
  201. 2160 for e=1 to n
  202. 2170 ap=sqr(ar(d,e)*ar(d,e)+ai(d,e)*ai(d,e))
  203. 2180 if ap=0 then 2330
  204. 2190 if ar(d,e)<>0 then 2230
  205. 2200 if ai(d,e)<0 then pa=-(NULL)/2
  206. 2210 if ai(d,e)>0 then pa=(NULL)/2
  207. 2220 goto 2320
  208. 2230 if ai(d,e)<>0 then 2270
  209. 2240 if ar(d,e)<0 the npa=(NULL)
  210. 2250 if ar(d,e)>0 thenpa=0
  211. 2260 goto 2320
  212. 2270 pa=atn(ai(d,e)/ar(d,e))
  213. 2280 if ar(d,e)>0 and ai(d,e)>0 thenpa=pa
  214. 2290 if ar(d,e)<0 and ai(d,e)>0 thenpa=(NULL)+pa
  215. 2300 if ar(d,e)<0 and ai(d,e)<0 thenpa=pa-(NULL)
  216. 2310 if ar(d,e)>0 and ai(d,e)<0 thenpa=pa
  217. 2320 ar(d,e)=ap: ai(d,e)=pa
  218. 2330 next e,d
  219. 2350 return
  220. 2360 :
  221. 2370 rem amplitude normalization routine
  222. 2380 rem find largest amplitude
  223. 2390 la=ar(1,1)
  224. 2400 for d=1 to n
  225. 2410 for e=1 to n
  226. 2420 ifar(d,e)>lathenla=ar(d,e)
  227. 2430 next e,d
  228. 2450 if ip$="r" then 2820
  229. 2460 print"the largest amplitude is";la
  230. 2470 rem check amplitude distribution
  231. 2480 p8=0: p6=0: p4=0: p2=0: p1=0
  232. 2490 n8=.8*la:n6=.6*la:n4=.4*la:n2=.2*la:n1=.1*la
  233. 2500 print" "
  234. 2510 for d=1 to n
  235. 2520 for e=1 to n
  236. 2530 if ar(d,e)<=n8 then p8=p8+1
  237. 2540 if ar(d,e)<=n6 then p6=p6+1
  238. 2550 if ar(d,e)<=n4 then p4=p4+1
  239. 2560 if ar(d,e)<=n2 then p2=p2+1
  240. 2570 if ar(d,e)<=n1 then p1=p1+1
  241. 2580 next e,d
  242. 2600 p8=int(p8/1024*100):p6=int(p6/1024*100):p4=int(p4/1024*100)
  243. 2610 p2=int(p2/1024*100):p1=int(p1/1024*100)
  244. 2620 print"80% of the max. amp is";n8
  245. 2630 print p8;"% of amp's are equal or smaller"
  246. 2640 print" "
  247. 2650 print"60% of the max. amp is";n6
  248. 2660 printp6;"% of amp's are equal or smaller"
  249. 2670 print" "
  250. 2680 print"40% of the max. amp is";n4
  251. 2690 printp4;"% of amp's are equal or smaller"
  252. 2700 print" "
  253. 2710 print"20% of the max. amp is";n2
  254. 2720 printp2;"% of amp's are equal or smaller"
  255. 2730 print" "
  256. 2740 print"10% of the max. amp is";n1
  257. 2750 printp1;"% of amp's are equal or smaller"
  258. 2760 print" "
  259. 2770 print"time is  ";ti$
  260. 2780 print" "
  261. 2790 input"what no. do you want to normalize to";nl
  262. 2800 print" "
  263. 2810 input"which means what % of data truncated";pt
  264. 2820 ifip$="r"thennl=la
  265. 2830 for d=1 to n
  266. 2840 for e=1 to n
  267. 2850 ar(d,e)=ar(d,e)/nl: if ar(d,e)>1 then ar(d,e)=1
  268. 2860 next e,d
  269. 2880 return
  270. 2890 :
  271. 2900 rem plot routine
  272. 2910 rem read in plot characters
  273. 2920 if hp$="n" then 3720
  274. 2930 print"...defining plot characters"
  275. 2940 for i=1 to 8:e(i)=2^i-1:e(i+8)=e(i)*2^(8-i):next
  276. 2950 for i=1 to 11: for j=1 to 15: for k=1 to 8
  277. 2960 ch=0: if k<(13-i) and k>(8-i) then ch=e(j)
  278. 2970 ch$(i,j)=ch$(i,j)+chr$(ch)
  279. 2980 next k,j,i
  280. 3710 :
  281. 3720 open 2,4,5
  282. 3730 open 6,4,6:print#6,chr$(21)
  283. 3740 rem main plotting loop
  284. 3750 tb=0
  285. 3760 nf=2: if hp$="n" then nf=1
  286. 3770 for d=1 to n
  287. 3780 for f=1 to nf
  288. 3790 for e=1 to n
  289. 3800 rem determine which character
  290. 3810 rem find character height
  291. 3820 ht=8*ar(d,e)
  292. 3830 dh=ht-int(ht)
  293. 3840 if dh<.125/2 then ht=int(ht):goto3860
  294. 3850 ht=int(ht)+1
  295. 3860 if f=nf then ar(d,e)=ht/8*nl
  296. 3870 rem find phase
  297. 3880 cp=pq*ai(d,e)/(2*(NULL))
  298. 3890 qp=int(pq/2)
  299. 3900 ifabs(cp)-int(abs(cp))>.5then3950
  300. 3910 if cp<0 then cp=int(cp):if f=nf then ai(d,e)=(cp+1)*(NULL)/qp
  301. 3920 if cp>=0 then cp=int(cp)+1:if f=nf then ai(d,e)=(cp-1)*(NULL)/qp
  302. 3930 if hp$="n"then 4290
  303. 3940 goto 3980
  304. 3950 if cp<0 then cp=int(cp)-1:if f=nf then ai(d,e)=(cp+1)*(NULL)/qp
  305. 3960 if cp>0then cp=int(cp)+2:if f=nfthen ai(d,e)=(cp-1)*(NULL)/qp
  306. 3970 if hp$="n" then 4290
  307. 3980 if f<>1 then 4230
  308. 3990 if ht=0 then print#4,"  ";chr$(141);:goto4260
  309. 4000 if abs(cp)>2then 4130
  310. 4010 if cp>-2 then 4050
  311. 4020 print#2,ch$(3,ht):print#4,chr$(254);chr$(141);
  312. 4030 tb=tb+1:print#2,ch$(11,ht):print#4,tab(tb);chr$(254);chr$(141);
  313. 4040 goto 4260
  314. 4050 if cp>1then 4090
  315. 4060 print#2,ch$(2,ht):print#4,chr$(254);chr$(141);
  316. 4070 tb=tb+1:print#2,ch$(10,ht):print#4,tab(tb);chr$(254);chr$(141);
  317. 4080 goto 4260
  318. 4090 print#2,ch$(1,ht):print#4,chr$(254);chr$(141);
  319. 4100 tb=tb+1:print#2,ch$(9,ht):print#4,tab(tb);chr$(254);chr$(141);
  320. 4110 tb=tb+1:print#2,ch$(9,ht):print#4,tab(tb);chr$(254);chr$(141);
  321. 4120 goto 4260
  322. 4130 if cp<0 then 4180
  323. 4140 cp=cp+11-cp*2
  324. 4150 print#4," ";chr$(141);
  325. 4160 tb=tb+1:print#2,ch$(cp,ht):print#4,tab(tb);chr$(254);chr$(141);
  326. 4170 goto4260
  327. 4180 if cp>0 then 4230
  328. 4190 cp=abs(cp)+1
  329. 4200 print#2,ch$(cp,ht):print#4,chr$(254);chr$(141);
  330. 4210 tb=tb+1:print#4,tab(tb);" ";chr$(141);
  331. 4220 goto 4260
  332. 4230 if ht=0 then 3990
  333. 4240 if ht<8 then ht=ht+8
  334. 4250 goto 3990
  335. 4260 if e=n then print#4,chr$(141);chr$(13);
  336. 4270 tb=2*e:if e=n then tb=0
  337. 4280 if e<>n then print#4,tab(tb);
  338. 4290 next e,f,d
  339. 4320 return
  340. 4330 :
  341. 4340 rem find real and imaginary components from truncated & quantized data
  342. 4350 for d=1 to n: for e=1 to n
  343. 4370 ap=ar(d,e)
  344. 4380 ar(d,e)=ap*cos(ai(d,e))
  345. 4390 ai(d,e)=ap*sin(ai(d,e))
  346. 4400 next e,d
  347. 4420 return
  348. 4430 :
  349. 4440 rem gray scale picture routine
  350. 4450 print#4," "
  351. 4460 for d=1 to n
  352. 4470 print#4,tab(23);
  353. 4480 for e=1 to n
  354. 4490 in=sqr(ar(d,e)*ar(d,e)+ai(d,e)*ai(d,e))
  355. 4500 af=1/sqr(ir*ir+ii*ii):in=af*in
  356. 4510 if in>1 then in=1
  357. 4520 if in<.1 thenprint#4," ";:goto4590
  358. 4530 if in<.2 thenprint#4,"'";:goto4590
  359. 4540 if in<.4 thenprint#4,":";:goto4590
  360. 4550 if in<.6 thenprint#4,"*";:goto4590
  361. 4560 if in<.8 thenprint#4,"[166]";:goto4590
  362. 4570 if in<.9 thenprint#4,"*[146]";:goto4590
  363. 4580 if in<=1 thenprint#4," [146]";
  364. 4590 if e=n then print#4,chr$(13);
  365. 4600 next e,d
  366. 4620 print#4," "
  367. 4630 print#4,"input real amplitude:";ir
  368. 4640 print#4,"input imaginary amplitude:";ii
  369. 4650 print#4,"phase angle:";ip$
  370. 4660 print#4,"no. of phase quantization levels:";pq
  371. 4670 if ip$="r" then 4700
  372. 4680 print#4,"transform normalized with...";nl
  373. 4690 print#4,"percent of transformed amp's truncated...";pt
  374. 4700 return
  375.