REM Cylindrical anamorphosis
REM Paste the transformed picture on the lateral fase of a cylinder.
OPTION ARITHMETIC NATIVE
SET COLOR MODE "NATIVE"        
GLOAD "zenkouji.JPG"              
ASK PIXEL SIZE (0,0; 1,1) a,b 
DIM c(a,b)                     
ASK PIXEL ARRAY (0,1) c                  ! 
CLEAR    
LET h = b * 5                            ! Viewpoint height
LET d = h * 1.0                          ! Horizontal distance between the cylinder and the viewpoint
LET m = 0                                ! Bottom margin
LET e=(b/a)*(d/h)*1.00       ! The ratio of the position of the point that is currectly reflected on the virtual surface.
LET k=1/(2+1/e)              ! The ratio of the point where the vitual surface is gererated.
LET r0 =SQR(1/(2*k-k^2))*(a/2)*1.00      ! The radius of the cylinder 
LET width=1000                           ! Width of the picture to be obtained.
LET HEIGHT=800                           ! Height of the picture to be obtained
LET top=200                              ! The position where the center of a cylinder to be put.
SET BITMAP SIZE width+1,height+1 
LET bottom=-height+top    
SET WINDOW -width/2, width/2, bottom , top   
DRAW grid(100,100)  
DRAW circle WITH SCALE(r0)   ! ‰~“›
LET d2=d-r0+r0*k                             ! Dsiatnce to the vitual surface 
SET POINT STYLE 1   
FOR x=-width/2 TO width/2 
   LET phi=0
   LET psi=0
   FOR y=bottom TO top 
      LET r=SQR(x^2+y^2)                     ! The absolute of the point to be painted
      IF r>r0 THEN
         LET t=ANGLE(x,y) +PI/2                
         IF t>PI THEN LET t=t-2*PI           ! The argument of the point to be painted
         FOR i=1 TO 2
            LET alpha=(t+psi+phi)/2  
            LET phi=ASIN(r0*SIN(alpha)/d) 
            LET psi=ASIN(r0*SIN(alpha)/r)
         NEXT i
         IF ABS(alpha)<=PI/2 THEN
            WHEN EXCEPTION IN 
               LET rr=r*SIN(alpha-psi)/SIN(alpha)
               LET dd=d*SIN(alpha-phi)/SIN(alpha)
            USE
               LET rr=r-r0
               LET dd=d-r0
            END WHEN
            LET y0=h/(1+dd/rr)
            LET x0=d2*TAN(phi)                       ! The x-coordinate of the original picuture
            LET y0=h-d2*(h-y0)/(d-r0*COS(alpha-phi)) ! The y-coordinate of the original picuture
            LET i=ROUND(x0+ (a+1)/2,0)               ! The index corresponding to x0
            LET j=b-ROUND(y0-m)                      ! The index corresponding to y0
            IF i>0 AND i<=a AND j>0 AND j<=b THEN
               SET POINT COLOR c(i,j)
               PLOT POINTS : x,y  
            END IF
         END if
      END if
   NEXT y
NEXT x
END
