!This code just for study, It doesn't work without certain input files.

finish
/clear,nostart

/INPUT,'Analysis','inp','','',0
!----------------normalized  curve start--------------------------------
 fix_p2z=1    ! 1: fix p2z=p3z, 0: free p2z
 f_line=end_line-start_line+1
 *dim,tbzr,table,f_line,1
 *dim,zrt,array,f_line,3

 *tread,tbzr(1,0),'%Filename1%',inp,,start_line-1
 *voper,zrt(1,1),1,mult,tbzr(1,0)
 *voper,zrt(1,2),1,mult,tbzr(1,1)
  cpts=4
 *dim,cpzr,array,cpts,2
!calculate the sum distance
 dist=0
 *do,i,1,f_line
 j=i+1
 dist_z=zrt(j,1)-zrt(i,1)
 dist_r=zrt(j,2)-zrt(i,2)
 dist=SQRT(dist_z**2+dist_r**2)
 sum=sum+dist
 *if,j,eq,f_line,exit
 *enddo

!calculate the normalized distance ratio
 dist=0
 zrt(1,3)=0
 zrt(f_line,3)=1
 *do,i,2,f_line-1
 j=i-1
 dist_z=zrt(j,1)-zrt(i,1)
 dist_r=zrt(j,2)-zrt(i,2)
 dist=dist+SQRT(dist_z**2+dist_r**2)
 zrt(i,3)=dist/sum
 *enddo
!----------------normalized curve end--------------------------------
!----------------fit curve start--------------------------------
/inquire,fil_ex,exist,control,pnt
*if,fil_ex,ne,1,then
!Using cubic bezier curve to fit the back-disk
!Create fuction: a1*P1+b1*P2+c1=0 and a2*P1+b2*P2+c2=0
!solve above fuction to get the control point P1 and P2
!Now we define the factor a, b, c
 cpzr(1,1)=zrt(1,1)
 cpzr(cpts,1)=zrt(f_line,1)
 cpzr(1,2)=zrt(1,2)
 cpzr(cpts,2)=zrt(f_line,2)
 *if,fix_p2z,eq,0,then    !free p2z
*set,a1,0
*set,a2,0
*set,b1,0
*set,b2,0
*set,c1,0
*set,c2,0
*set,p1z,0
*set,p2z,0
 *do,i,1,f_line
 t=zrt(i,3)
 p0z=cpzr(1,1)
 p3z=cpzr(cpts,1)
 zt=zrt(i,1)
 a1=a1+3*(1-t)**2*t*(3*(1-t)**2*t)
 b1=b1+3*(1-t)*t**2*(3*(1-t)**2*t)
 c1=c1+((1-t)**3*p0z+t**3*p3z-zt)*(3*(1-t)**2*t)
 a2=a2+3*(1-t)**2*t*(3*(1-t)*t**2)
 b2=b2+3*(1-t)*t**2*(3*(1-t)*t**2) 
 c2=c2+((1-t)**3*p0z+t**3*p3z-zt)*(3*(1-t)*t**2)
 *enddo
 P1z=(b1*c2-c1*b2)/(a1*b2-b1*a2)
 P2z=(a1*c2-c1*a2)/(b1*a2-a1*b2)
 cpzr(2,1)=P1z
 cpzr(3,1)=P2z
*set,a1,0
*set,a2,0
*set,b1,0
*set,b2,0
*set,c1,0
*set,c2,0
*set,p1r,0
*set,p2r,0
 *do,i,1,f_line
 t=zrt(i,3)
 p0r=cpzr(1,2)
 p3r=cpzr(cpts,2)
 rt=zrt(i,2)
 a1=a1+3*(1-t)**2*t*(3*(1-t)**2*t)
 b1=b1+3*(1-t)*t**2*(3*(1-t)**2*t)
 c1=c1+((1-t)**3*p0r+t**3*p3r-rt)*(3*(1-t)**2*t)
 a2=a2+3*(1-t)**2*t*(3*(1-t)*t**2)
 b2=b2+3*(1-t)*t**2*(3*(1-t)*t**2) 
 c2=c2+((1-t)**3*p0r+t**3*p3r-rt)*(3*(1-t)*t**2)
 *enddo
 P1r=(b1*c2-c1*b2)/(a1*b2-b1*a2)
 P2r=(a1*c2-c1*a2)/(b1*a2-a1*b2)
 cpzr(2,2)=P1r
 cpzr(3,2)=P2r

 *else        !fix p2z
 *set,a1,0
 *set,c1,0
 *set,p1z,0
 p0z=cpzr(1,1)
 p2z=cpzr(cpts,1)
 p3z=cpzr(cpts,1)
 *do,i,1,f_line
 t=zrt(i,3)
 zt=zrt(i,1)
 a1=a1+3*(1-t)**2*t*(3*(1-t)**2*t)
 par1=(1-t)**3*p0z+3*(1-t)*t**2*p2z
 c1=c1+(par1+t**3*p3z-zt)*(3*(1-t)**2*t)
 *enddo
 P1z=(-1)*c1/a1
 cpzr(2,1)=P1z
 cpzr(3,1)=P2z
 *set,a1,0
 *set,a2,0
 *set,b1,0
 *set,b2,0
 *set,c1,0
 *set,c2,0
 *set,p1r,0
 *set,p2r,0
 *do,i,1,f_line
 t=zrt(i,3)
 p0r=cpzr(1,2)
 p3r=cpzr(cpts,2)
 rt=zrt(i,2)
 a1=a1+3*(1-t)**2*t*(3*(1-t)**2*t)
 b1=b1+3*(1-t)*t**2*(3*(1-t)**2*t)
 c1=c1+((1-t)**3*p0r+t**3*p3r-rt)*(3*(1-t)**2*t)
 a2=a2+3*(1-t)**2*t*(3*(1-t)*t**2)
 b2=b2+3*(1-t)*t**2*(3*(1-t)*t**2) 
 c2=c2+((1-t)**3*p0r+t**3*p3r-rt)*(3*(1-t)*t**2)
 *enddo
 P1r=(b1*c2-c1*b2)/(a1*b2-b1*a2)
 P2r=(a1*c2-c1*a2)/(b1*a2-a1*b2)
 cpzr(2,2)=P1r
 cpzr(3,2)=P2r
 *endif  !end fix_p2z

*cfopen,control,pnt
*vwrite,'Z      R'
%C
*vwrite,cpzr(1,1),cpzr(1,2)
(3x,f11.8,4x,f11.8)
*cfclos
*endif     !end fitting