Nose Shapes with Bluntness

December 23, 2010

This post shape formulation for three nose shapes:

  • Tangent ogive
  • Conic ogive
  • Double arc ogive
  • Tangent ogive

    dimensions of Tangent ogive for analysis

    \displaystyle \begin{array}{c}b+\left(R-b\right)\sin\theta=L_{n} \\ b\cos\theta+\left(R-R\cos\theta\right)=R_{n}\end{array}

    solving

    \displaystyle \begin{array}{c}R=\frac{1}{2}\frac{L_{n}^{2}-2bL_{n}+R_{n}^{2}}{R_{n}-b}\\ \theta=\arctan\left(\frac{2\left(L_{n}-b\right)\left(R_{n}-b\right)}{\left(R_{n}-L_{n}\right)\left(2b-\left(L_{n}+R_{n}\right)\right)}\right)\end{array}

    \left(Center,Radius,\left[\theta_{1},\theta_{2}\right]\right) descriptions of two arcs are as follows:

    \displaystyle Arc_{1}\rightarrow\left(b,0\right),b,\left[\frac{\pi}{2}+\theta,\pi\right]\\Arc_{2}\rightarrow\left(L_{n},R_{n}-R\right),R,\left[\frac{\pi}{2},\frac{\pi}{2}+\theta\right]

    where R is the radius of nose, b is the radius of bluntness, and \theta is arc-angle of bluntness

    A plot of family of blunt tangent ogives for varying bluntness radius are given below

    a family of tangent ogives with varying bluntness

    PDF animation is generated in Asymptote.

    Tangent ogive with varying nose bluntness

    Conic Ogive

    Conic ogive dimensions for analysis

    For conic nose there is only one parameter (\theta) to be determined. Thus only one equation is needed. Writing y equation

    \left(L_{n}-b\right)\sin\theta=R_{n}\cos\theta-b\\R_{n}\cos\theta+\left(b-L_{n}\right)\sin\theta=b

    solving for \theta

    \displaystyle e^{i\theta}=\frac{bR_{n}-K\left(b-L_{n}\right)+i\left(b\left(b-L_{n}\right)+KR_{n}\right)}{R_{n}^{2}+\left(b-L_{n}\right)^{2}}

    family of conic ogive for varying bluntness

    PDF animation is generated by Asymptote

    Conic ogive with bluntness

    Double arc Ogive

    Double arc nose dimensions for analysis

    Direct solution of x-y equations for double-arc nose

    \displaystyle \overline{BA}+\overline{BG}=R_{1}-\sqrt{R_{1}^{2}-H^{2}}+x=L_{n}\\\overline{BD}+\overline{HF}=H+R_{2}-\sqrt{R_{2}^{2}-x^{2}}=R_{n}

    results in complicated expressions, so using

    \displaystyle \overline{BC}=R_{1}\sin2\theta\\\overline{BA}=R_{1}-\overline{BC}=R_{1}-R_{1}\sin2\theta\\\overline{BG}=R_{2}\sin2\theta
    and
    \displaystyle \overline{BD}=H=R_{1}\cos2\theta\\\overline{HF}=h=R_{2}-R_{2}\cos2\theta

    results in x-y equations

    \displaystyle R_{1}-R_{1}\sin2\theta+R_{2}\sin2\theta=L_{n}\\R_{1}\cos2\theta+R_{2}-R_{2}\cos2\theta=R_{n}

    simplifying these equations

    \displaystyle \sin2\theta\left[R_{2}-R_{1}\right]=L_{n}-R_{1}\\\cos2\theta\left[R_{1}-R_{2}\right]=R_{n}-R_{2}

    summing squares

    \displaystyle 2L_{n}R_{1}+2R_{n}R_{2}-2R_{2}R_{1}=L_{n}^{2}+R_{n}^{2}

    if ratio of radii information \lambda=\frac{R_{1}}{R_{2}} is imposed

    \displaystyle 0=2\lambda R_{2}^{2}-2\left(L_{n}\lambda+R_{n}\right)R_{2}+L_{n}^{2}+R_{n}^{2}

    if nose fineness ratio f=\frac{L_{n}}{R_{n}} parameter is used

    \displaystyle 0=2\lambda R_{2}^{2}-2R_{n}\left(f\lambda+1\right)R_{2}+R_{n}^{2}\left(1+f^{2}\right)

    After determining R_{2} from this equation, R_{1} is computed via R_{1}=\lambda R_{2}. To compute \theta

    \displaystyle \sin2\theta=\frac{L_{n}-R_{1}}{R_{2}-R_{1}}\\\cos2\theta=\frac{R_{2}-R_{n}}{R_{2}-R_{1}}\\e^{i2\theta}=\frac{\left(R_{2}-R_{n}\right)+i\left(L_{n}-R_{1}\right)}{R_{2}-R_{1}}

    Using XFOIL and Automating via Python subprocess module

    December 21, 2010

    XFOIL is an airfoil profile analysis program developed by Mark Drela. Usage is through command prompt where user issues from a set of commands to load/modify airfoil profiles, perform aerodynamic analysis of airfoil and store results as graphics and in tabular text format.

    A short walkthrough of typical session (from the XFOIL documentation)

    A basic command set from XFOIL documentation “sessions.txt” is given below

    a session fragment from XFOIL documentation
    LOAD e387.dat
    GDES      (enter GDES menu)        |
    CADD      (add points at corners)  |  These commands are optional,
     (accept default input)   |  and are recommended only for
     (accept default input)   |  Eppler and Selig airfoils
     (accept default input)   |  to give smoother LE shapes
     (return to Top Level)    |
    
    PANEL    (regenerate paneling since)
                (better panel node spacing is needed)
    OPER     (enter OPER menu)
    ALFA 5.0
    
  • Command “LOAD ../runs/e387.dat” loads the airfoil profile from file e387.dat

    plot window after loading airfoil profile

  • command “GDES” switches the XFOIL to geometry design stage
  • command “CADD” adds corner points to satisfy user specified angle threshold. After this command four arguments are asked to the user. The session given above accepts defaults (by “enter”ing)
  • command “PANEL” regenerates vertices
  • command “OPER” switches the program to aerodynamic analysis mode
  • command “ALFA 5.0” sets the angle of attack to 5 degrees (after this command analysis for the given AoA is performed immediately and plot window gives pressure coefficient C_p plot

    Cp plot for AoA=5 deg

  • Automating the execution

    To automate execution of XFOIL python subprocess package is used. XFOIL is opened as a process and and usual XFOIL prompt commands are fed to input stream of the process. Commands issued are same as the command sequence from XFOIL sample session document.

    # -*- coding: utf-8 -*-
    # FILE: runXFOIL_alphaSweep.py
    
    import subprocess as sp
    #import os
    import shutil
    import sys
    import string
    
    """
    class Popen(args, bufsize=0, executable=None,
                stdin=None, stdout=None, stderr=None,
                preexec_fn=None, close_fds=False, shell=False,
                cwd=None, env=None, universal_newlines=False,
                startupinfo=None, creationflags=0):
    """
    
    ps = sp.Popen(['xfoil.exe'],
                  stdin=sp.PIPE,
                  stdout=None,
                  stderr=None)
    
    def issueCmd(cmd,echo=True):
        ps.stdin.write(cmd+'\n')
        if echo:
            print cmd
    
    method = 2
    
    if method==1:
        res = ps.communicate( string.join(["load ../runs/e387.dat",
                                           "oper",
                                           "alfa 0.0",
                                           "cpwr cp_a0.dat",
                                           "hard",
                                           "alfa 2.0",
                                           "cpwr cp_a2.dat",
                                           "hard",
                                           "alfa 4.0",
                                           "cpwr cp_a4.dat",
                                           "hard",
                                           "alfa 6.0",
                                           "cpwr cp_a6.dat",
                                           "hard"],'\n') )
        print res
    elif method==2:
        issueCmd('load ../runs/e387.dat\n')
    
        """
        GDES              (enter GDES menu)        |
        CADD              (add points at corners)  |  These commands are optional,
                  (accept default input)   |  and are recommended only for
                  (accept default input)   |  Eppler and Selig airfoils
                  (accept default input)   |  to give smoother LE shapes
                  (return to Top Level)    |
    
        PANEL           (regenerate paneling since better panel node spacing is needed)
        """
        issueCmd('GDES')  # enter GDES menu
        issueCmd('CADD')  # add points at corners
        issueCmd('')      # accept default input
        issueCmd('')      # accept default input
        issueCmd('')      # accept default input
        issueCmd('')      # accept default input
        issueCmd('PANEL') # regenerate paneling
        issueCmd('SAVE e387coords.dat') # save panel geometry
        issueCmd('OPER')
        count = 0
        for aoa in [-10+0.2*i for i in range(101)]:
            issueCmd('ALFA %7.4f' % (aoa,))
            issueCmd('CPWR cpprof%02d.dat' % (count,))
            issueCmd('HARD')
            #print "resp:",ps.stdout.read()
            print "renaming",count
            shutil.copyfile('plot.ps','plot%02d.ps' % (count,))
            count += 1
    
        resp = ps.stdout.read()
        print "resp:",resp
    else:
        pass
    

    In this script, i have tried two methods to interact with XFOIL

  • communicate” method of process object
  • write” method of input stream of process object
  • My questionable (since this is my initial experimentation with subprocess module of Python) conclusions are:

  • It seems that communicate requires both input and output streams when constructing the process object and is suitable for one time communication (since it waits for the termination of program).
  • On the other hand using input/output streams, commands to executed process can be supplied to successively, allowing the usage scenario shape adaptively
  • i experience hangup problems when tried to read from output stream of process object, so later i nullified this stream
  • Postprocessing Results

    # FILE: processXFOILresults.py
    
    import string
    import pylab as pl
    import glob
    
    from matplotlib import rc
    
    # ------------------------------------------------------
    """
    font.family        : serif
    font.serif         : Times, Palatino, New Century Schoolbook, Bookman, Computer Modern Roman
    font.sans-serif    : Helvetica, Avant Garde, Computer Modern Sans serif
    font.cursive       : Zapf Chancery
    font.monospace     : Courier, Computer Modern Typewriter
    
    text.usetex        : true
    """
    USETEX = True
    if USETEX:
        rc('font',**{'family':'sans-serif',
                     'sans-serif':['Computer Modern Roman'],
                     'serif':'Computer Modern Roman',
                     'monospace':'Computer Modern Typewriter'})
        rc('text', usetex=True)
    
    # ------------------------------------------------------
    def line2floats(line):
        try:
            parts = string.split( line )
            floats = map( float, parts)
        except:
            parts = line[0:10], line[10:]
            floats = map( float, parts)
    
        return floats
    
    def readTwoColumnFile(fname):
        fid = file(fname,'rt')
        lines = fid.readlines()
        lines = lines[1:] # remove header
    
        data = map(line2floats, lines)
        data = pl.array(data)
        fid.close()
    
        return data
    
    def bisector(a,b):
        " assumes a and b are unit length"
        return (a+b)/pl.norm(a+b)
    
    def postprocessCpData(data,geo,newxcount):
        x = data[:,0]
        y = geo[:,1]
        Cp = data[:,1]
    
        n = data.shape[0]
    
        # compute finite difference of x to classify points as upper and lower airfoil surface
        dx = pl.diff(x)
        dy = pl.diff(y)
        L = pl.sqrt(dx**2+dy**2)
        Tx = dx/L
        Ty = dy/L
        Nx = -Ty
        Ny = +Tx
        T = pl.array((Tx,Ty))
        T = T.transpose()
        N = pl.array((Nx,Ny))
        N = N.transpose()
    
        midx = (x[0:n-1]+x[1:n])/2.0
        midy = (y[0:n-1]+y[1:n])/2.0
        midcp = (Cp[0:n-1]+Cp[1:n])/2.0
    
        Tnode = pl.zeros( (n,2), pl.float64)
        Nnode = pl.zeros( (n,2), pl.float64)
        for i in range(1,n-1):
            Tnode[i,:] = bisector( T[i-1,:], T[i,:] )
            Nnode[i,:] = bisector( N[i-1,:], N[i,:] )
        Tnode[0,:] = bisector( T[0,:], T[-1,:] )
        Tnode[-1,:] = bisector( T[0,:], T[-1,:] )
        Nnode[0,:] = bisector( N[i-1,:], N[i,:] )
        Nnode[-1,:] = bisector( N[i-1,:], N[i,:] )
    
        # determine (safe) limits of x for interpolation
        xmin = min( min(x[dx<0]),min(x[dx>=0]) )
        xmax = max( max(x[dx<0]),max(x[dx>=0]) )
    
        # re-compute lower and upper Cp at new abscissae
        if ChebyshevSpacing:
            xnew = pl.linspace(pl.pi, 0, newxcount)
            xnew = (pl.cos(xnew)+1)/2.0
        else:
            xnew = pl.linspace(xmin, xmax, newxcount)
    
        newCpUpper = pl.interp(xnew, pl.flipud(x[dx<0]), pl.flipud(Cp[dx<0]))     newCpLower = pl.interp(xnew, x[dx>=0], Cp[dx>=0])
    
        return (x,y,Cp,L,T,N,midx,midy,midcp,Tnode,Nnode,xnew,newCpUpper,newCpLower)
    
    # ------------------------------------------------------
    
    datadir = './plots2'
    
    CpFiles = glob.glob(datadir+'/dat/cpprof*.dat')
    geofile = datadir+'/dat/e387coords.dat'
    
    geodata = readTwoColumnFile(geofile)
    
    count = 0
    maxcount = 5
    ChebyshevSpacing = True
    
    postprocessResult = file('post.dat','wt')
    postprocessResult.write( "%12s %12s %12s %12s %12s %12s\n" % ("AoA","FX","FY","MY","FD","FL") )
    
    ForceMoment = []
    
    for CpFile in CpFiles:
        alpha = -10+0.2*count
        print "reading Cp file" , CpFile
        data = readTwoColumnFile(CpFile)
    
        x, y, Cp, L, T, N, midx, midy, midcp, Tnode, Nnode, xnew, newCpUpper, newCpLower = postprocessCpData(data, geodata, 100)
    
        # panel Cp x and y components
        midcpx = midcp*N[:,0]
        midcpy = midcp*N[:,1]
    
        dx = pl.diff(x)
        dy = pl.diff(y)
    
        forceX = pl.trapz(midcpx*L)
        forceY = pl.trapz(midcpy*L)
        momentY = pl.trapz((midx-0.25)*midcpy*L)
        AoA = alpha*pl.pi/180.
        forceD = +forceX*pl.cos(AoA)+forceY*pl.sin(AoA)
        forceL = -forceX*pl.sin(AoA)+forceY*pl.cos(AoA)
        postprocessResult.write( "%12.6f %12.6f %12.6f %12.6f %12.6f %12.6f\n" % (AoA,forceX,forceY,momentY,forceD, forceL) )
    
        ForceMoment.append( [AoA,forceX,forceY,momentY,forceD,forceL])
    
        # plot upper and lower profile parts in different colors
        pl.figure()
        pl.plot( x[dx<0] ,data[dx<0] ,'r.-',lw=2)     pl.plot( x[dx>=0],data[dx>=0],'b.-',lw=2)
        pl.legend(['upper $C_p$','lower $C_p$'])
        pl.title(r'$C_p @ \alpha -10$')
    
        # plot upper and lower Cp in different color
        pl.figure()
        pl.plot( xnew,newCpUpper,'r.-',lw=2)
        pl.plot( xnew,newCpLower,'b.-',lw=2)
        pl.legend(['upper $C_p$','lower $C_p$'])
        pl.title(r'$C_p (\alpha=%6.3f)$' % (alpha,))
    
        # plot Cp difference between upper and lower profiles
        pl.figure()
        pl.plot( xnew,newCpUpper-newCpLower,'r.-',lw=2)
        pl.title(r'$\Delta C_p (\alpha=%6.3f)$' % (alpha,))
        pl.ylim([-5,14])
    
        pl.savefig(datadir+'/fig%02d.pdf' % (count,))
    
        count += 1
        if count>maxcount:
            break
    
    postprocessResult.close()
    
    ForceMoment = pl.array(ForceMoment)
    
    pl.figure()
    pl.plot( ForceMoment[:,0], ForceMoment[:,1], 'd-', ForceMoment[:,0], ForceMoment[:,2], 'o-', mec=None)
    pl.xlabel('$x$')
    pl.ylabel('$F$')
    pl.legend( ('$F_x$','$F_y$') )
    pl.grid()
    pl.savefig('forceXY.pdf')
    
    pl.figure()
    pl.plot( ForceMoment[:,0], ForceMoment[:,4], 'd-', ForceMoment[:,0], ForceMoment[:,5], 'o-', mec=None)
    pl.xlabel('$x$')
    pl.ylabel('$F$')
    pl.legend( ('$F_D$','$F_L$') )
    pl.grid()
    pl.savefig('forceDL.pdf')
    
    pl.figure()
    pl.plot( ForceMoment[:,0], ForceMoment[:,3], 'd-', mec=None)
    pl.xlabel('$x$')
    pl.ylabel('$M_y$')
    pl.grid()
    pl.savefig('momentY.pdf')
    

    To integrate C_p on discrete geometry, panel and node tangents are computed by

    \displaystyle \delta\mathbf{v}_{i}=\mathbf{v}_{i+1}-\mathbf{v}_{i}\\\mathbf{T}_{i}^{panel}=\frac{\delta\mathbf{v}_{i}}{\left\Vert \delta\mathbf{v}_{i}\right\Vert }\\\mathbf{N}_{i}^{panel}=\left[-T_{iy}^{panel},T_{ix}^{panel}\right]\\\mathbf{T}_{i}^{node}=\frac{\mathbf{T}_{i}^{panel}+\mathbf{T}_{i+1}^{panel}}{\left\Vert \mathbf{T}_{i}^{panel}+\mathbf{T}_{i+1}^{panel}\right\Vert }\\\mathbf{N}_{i}^{node}=\left[-T_{iy}^{node},T_{ix}^{node}\right]

    Total coefficient is computed by summing forces over all panels

    \displaystyle \mathbf{F}=\sum_{i}\left(C_{p}\right)_{i}L_{i}\mathbf{N}_{i}^{panel}

    Coefficients in airfoil axis is converted to flow axis via

    \displaystyle C_{D}=C_{X}\cos\left(\alpha\right)+C_{Y}\sin\left(\alpha\right)\\C_{L}=-C_{X}\sin\left(\alpha\right)+C_{Y}\cos\left(\alpha\right)

    graphic outputs of this code (axial-normal, drag lift and pitching moment coefficient vs. AoA) are given below

    Axial and normal aerodynamic force coefficients

    Drag and lift aerodynamic force coefficients

    Aerodynamic pitching Moment Coefficient

    to gain insight to the analysis routine def postprocessCpData(data,geo,newxcount) shapes (~sizes) of constructed result arrays are given below. Note that panel quantities are one less than the those of node.

    >>> x.shape
    (160,)
    >>> y.shape
    (160,)
    >>> L.shape
    (159,)
    >>> T.shape
    (159, 2)
    >>> N.shape
    (159, 2)
    >>> midx.shape
    (159,)
    >>> midy.shape
    (159,)
    >>> midcpx.shape
    (159,)
    >>> midcpy.shape
    (159,)
    >>> Tnode.shape
    (160, 2)
    >>> Nnode.shape
    (160, 2)
    >>> xnew.shape
    (100,)
    >>> newCpUpper.shape
    (100,)
    >>> newCpLower.shape
    (100,)
    >>>
    

    Rendering and Animating Selig Airfoils in Asymptote

    December 12, 2010

    The UIUC Airfoil Data Site contains a collection of airfoil profiles. Profiles are listed online and whole collection is downloadable (take coord_seligFmt.zip in page or direct link coord_seligFmt.zip)

    Airfoil profiles are contained in individual text files with one descriptive header line and lines of (x,y) coordinate pairs. To exemplify AG38.dat contents are partially listed below

    AG38
         0.999999    0.004704
         0.993621    0.005696
         0.983682    0.007241
         0.972066    0.009047
         0.959633    0.010981
         ...
         ...
         0.994806    0.000000
         1.000000    0.000001
    

    Coordinates start from aft of airfoil to fore of airfoil by upper part and then goes back to aft following lower part.

    Asymptote is a vector rendering library akin to MPOST. The Asymptote source code given below renders the airfoil profile in file “ultimate.dat”.

    size(10cm,0);
    
    write("started");
    
    file fin=input("./SeligAirfoils/ultimate.dat");
    
    string header = fin;
    write("header is ", header);
    
    real coords[][];
    
    guide allprof, upprof,loprof;
    
    int count = 0;
    while(true){
    	tmp = fin.line();
    	coords.push( tmp );
    	allprof = allprof .. (tmp[0],tmp[1]);
    	write("coords :", coords[count]);
    	++count;
    	if(eof(fin))
    		break;
    }
    
    draw(allprof .. cycle, black+1bp);
    
    write("coords length : ", coords.length);
    write("ended");
    close(fin);
    
  • To open a file for reading, command “input” is used (line 5)
  • one line of string is read by assignment “string=file” as described in Asymptote manual “Strings are also read by assignment, by reading characters up to but not including a newline“. Header text is read this way (line 7)
  • array of arrays to hold x,y coordinates is declared in line 10
  • all coordinates are read in while loop (lines 15-23) and appended to both “coords” array and “Asymtote guide” (akin to “path”) “allprof”.
  • finally path is rendered in line 25
  • Output of this code is readAndPlotSeligProfile

    Simplest render of ultimate airfoil from Selig Airfoil Database

    To split upper and lower parts of profile and determine camber line from those two, above code is modified as follows

    size(20cm,0);
    unitsize(1cm);
    
    write("started");
    
    string profname = "./SeligAirfoils/naca643218.dat";
    
    file fin=input(profname);
    
    string header = fin;
    write("header is ", header);
    
    guide allprof, upprof,loprof;
    
    real coord[];
    
    coord = fin.line();
    allprof = (coord[0],coord[1]);
    upprof = (coord[0],coord[1]);
    
    real lastx, currentx;
    
    lastx = coord[0];
    
    int count = 0;
    while(true){
    	coord = fin.line();
    	allprof = allprof .. (coord[0],coord[1]);
    	currentx = coord[0];
    	if(lastx-currentx>0)
    		upprof = upprof .. (coord[0],coord[1]);
    	else
    		loprof = loprof .. (coord[0],coord[1]);
    	lastx = currentx;
    	if(eof(fin)) break;
    }
    
    /* align top and bottom profiles in same direction ~ from rear to front */
    guide upprofp = upprof;
    guide loprofp = reverse(loprof);
    
    int i;
    guide chordline;
    for(i=0;i
    	chordline = chordline .. (point(upprofp,i)+point(loprofp,i))/2;
    }
    
    //draw(allprof .. cycle, black+1bp);
    draw(upprof , red+3bp+opacity(0.5), ArcArrow( DefaultHead, size=8, filltype=NoFill ));
    draw(loprof , green+3bp+opacity(0.5), ArcArrow( DefaultHead, size=8, filltype=NoFill ));
    
    draw( chordline, gray+1bp+dashed);
    
    real flapfrac = 0.3;
    
    /* determine the curve parameter for the given arclength fraction "flapfrac" */
    real uptime = reltime(upprofp, flapfrac);
    real lotime = reltime(loprofp, flapfrac);
    
    /* determine junction point position and tangent */
    pair juncposup = point(upprofp, uptime);
    pair junctanup = dir(upprofp, uptime);
    
    pair juncposlo = point(loprofp, lotime);
    pair junctanlo = dir(loprofp, lotime);
    
    /* split upper profile to front and rear parts */
    guide upproffront = subpath( upprofp, uptime, length(upprofp));
    guide upprofrear = subpath( upprofp, 0, uptime);
    
    /* split lower profile to front and rear parts */
    guide loproffront = subpath( loprofp, lotime, length(loprofp));
    guide loprofrear = subpath( loprofp, 0, lotime);
    
    /* determine rotation center */
    path auxpathup = juncposup -- juncposup + 99 * I * junctanup;
    path auxpathlo = juncposlo -- juncposlo - 99 * I * junctanlo;
    
    pair rotcenter;
    rotcenter = (juncposlo+juncposup)/2;
    //rotcenter = intersectionpoint( auxpathup, auxpathlo );
    
    transform T = shift(0,-0.3);
    
    label( header, (0.5,-0.5));
    
    draw( T*upproffront, red+dashed+1bp);
    draw( T*upprofrear, red+white+dashdotted+3bp);
    
    draw( T*loproffront, green+dashed+1bp);
    draw( T*loprofrear, green+white+dashdotted+3bp);
    
    draw( juncposlo -- rotcenter -- juncposup, blue);
    
    dot( Label("C"), rotcenter, magenta+3bp);
    
    write("ended");
    close(fin);
    
  • this code renders the profile “naca643218.dat” (line 6)
  • first coordinate is read and upper/all profiles are initialized by first coordinate read (line 17-19)
  • while loop is modified so as to notice whether x coordinate is increasing or decreasing, upon which coordinates read in while loop are appended to upper or lower profiles (lines 21-36)
  • align top and bottom profiles in same direction ~ from rear to front (lines 39-40)
  • chordline is constructed by traversing the all points of lower and upper profiles and computing the mean of two (lines 39-40). Note that this assumes that both profiles have equal or enough number of vertices
  • lower and upper profiles are rendered with red and green pens (lines 49-50)
  • curve parameter for given flap position “flapfrac” is computed (lines 57-58)
  • flap junction positions and tangents are computed (lines 61-65)
  • upper and lower profiles are partitioned to front (before flap position) and rear parts (after flap position) (lines 68-73)
  • rotation center for airfoil is computed by two methods (one is commented) (lines 75-81)
  • partioned airfoil profile is rendered below the first rendering (lines 83-95)
  • This code results in the readProfileSimple_save1 and converted to PNG given below

    Asymptote rendering of lower and upper parts of naca643218 from Selig Airfoil Database

     

    Animation of Flap Motion

    Previous Asymptote code is modified to animate flap motion. Flap is rotated about the rotation center

    import animate;
    
    settings.tex="pdflatex";
    settings.keep = true;
    
    //size(20cm,0);
    unitsize(10cm);
    
    write("started");
    
    string profname = "../SeligAirfoils/naca643218.dat";
    
    file fin=input(profname);
    
    string header = fin;
    write("header is ", header);
    
    guide allprof, upprof,loprof;
    
    real coord[];
    
    coord = fin.line();
    allprof = (coord[0],coord[1]);
    upprof = (coord[0],coord[1]);
    
    real lastx, currentx;
    
    lastx = coord[0];
    
    int count = 0;
    while(true){
    	coord = fin.line();
    	allprof = allprof .. (coord[0],coord[1]);
    	currentx = coord[0];
    	if(lastx-currentx>0)
    		upprof = upprof .. (coord[0],coord[1]);
    	else
    		loprof = loprof .. (coord[0],coord[1]);
    	lastx = currentx;
    	if(eof(fin)) break;
    }
    
    // insert last point of upper profile to beginning of lower profile
    loprof = point(upprof, length(upprof)) .. loprof;
    
    /* align top and bottom profiles in same direction ~ from rear to front */
    guide upprofp = upprof;
    guide loprofp = reverse(loprof);
    
    int i;
    guide chordline;
    for(i=0;i
    	chordline = chordline .. (point(upprofp,i)+point(loprofp,i))/2;
    }
    
    real flapfrac = 0.3;
    real shorten = 0.01;
    
    /* determine the curve parameter for the given arclength fraction "flapfrac" */
    real uptime = reltime(upprofp, flapfrac);
    real uptime_pos = reltime(upprofp, flapfrac+shorten);
    real uptime_neg = reltime(upprofp, flapfrac-shorten);
    
    real lotime = reltime(loprofp, flapfrac);
    real lotime_pos = reltime(loprofp, flapfrac+shorten);
    real lotime_neg = reltime(loprofp, flapfrac-shorten);
    
    /* determine junction point position and tangent */
    pair juncposup = point(upprofp, uptime);
    pair junctanup = dir(upprofp, uptime);
    
    pair juncposlo = point(loprofp, lotime);
    pair junctanlo = dir(loprofp, lotime);
    
    /* split upper profile to front and rear parts */
    guide upproffront = subpath( upprofp, uptime_pos, length(upprofp));
    guide upprofrear = subpath( upprofp, 0, uptime_neg);
    
    /* split lower profile to front and rear parts */
    guide loproffront = subpath( loprofp, lotime_pos, length(loprofp));
    guide loprofrear = subpath( loprofp, 0, lotime_neg);
    
    /* determine rotation center */
    path auxpathup = juncposup -- juncposup + 99 * I * junctanup;
    path auxpathlo = juncposlo -- juncposlo - 99 * I * junctanlo;
    
    pair rotcenter;
    rotcenter = (juncposlo+juncposup)/2;
    //rotcenter = intersectionpoint( auxpathup, auxpathlo );
    
    transform T = shift(0,-0.3);
    
    animation a;
    
    pen airfoilpen = blue+2bp;
    
    for(i=-10;i<10;++i){
    	write("doing step ",i);
    
    	T = rotate( i, rotcenter );
    
    	draw( upproffront, airfoilpen);
    	draw( T*upprofrear, airfoilpen);
    
    	draw( loproffront, airfoilpen);
    	draw( T*loprofrear, airfoilpen);
    
    	/* at junction point */
    	// draw( juncposup{junctanup*I} .. {junctanlo*I}juncposlo, blue );
    
    	/* at shortened junction points */
    	if(false){
    		// S-shaped
    		draw( point(loproffront,0){dir(loproffront,0)} .. {dir(upproffront,0)}point(upproffront,0), airfoilpen );
    		draw( T*(point(loprofrear,length(loprofrear)) {dir(loprofrear,length(loprofrear))} .. {dir(upprofrear,length(upprofrear))} point(upprofrear,length(upprofrear))), airfoilpen );
    	} else {
    		// (-shaped
    		draw( point(loproffront,0) {-I*dir(loproffront,0)} .. {-I*dir(upproffront,0)} point(upproffront,0), airfoilpen );
    		draw( T*(point(loprofrear,length(loprofrear)) {-I*dir(loprofrear,length(loprofrear))} .. {-I*dir(upprofrear,length(upprofrear))} point(upprofrear,length(upprofrear))), airfoilpen );
    	}
    
    	/* junction points connected */
    	//draw( juncposlo -- rotcenter -- juncposup, blue);
    
    	draw( chordline, gray+1bp+dashed);
    
    	dot( Label("C"), rotcenter, magenta+3bp);
    
    	a.add( currentpicture );
    
    	erase( currentpicture );
    }
    
    label(a.pdf("controls",multipage=false));
    
    write("ended");
    close(fin);
    

    this file results in the readProfileSimple_animate_forWordpress_save1 from which the following animation is generated.

    To generate animation each frame of PDF is stored separately or together (via settings.keep=true or multipage=true) and then converted to animation using ImageMagick. _readProfileSimple_animate_forWordpress_multipage is generated via “label(a.pdf(“controls”,multipage=true));”

    multipage PDF is converted to GIF animation via following ImageMagick commands

    rem makeAnim.bat
    rem !! make double percents %% single percent % when directly issuing commands
    rem generate PNG for each page of multipage PDF
    convert -filter Mitchell -resize 400x400 -density 300x300 _readProfileSimple_animate_forWordpress_multipage.pdf _readProfileSimple_animate_forWordpress_multipage%%02d.png
    
    rem convert/combine generated PNG to animated GIF
    convert _readProfileSimple_animate_forWordpress_multipage%%02d.png[0-19] -loop 0 airfoilgifanim.gif
    

    Selif Airfoil Flap Animation

    Statistical and Geometric Interpretations of Kalman Filter

    September 5, 2010

    In this post, geometrical interpretation(s) of Kalman filter is exemplified together with some visualizations.

    Bayesian View

    If current state (prior) estimate is described by Gaussian PDF

    \displaystyle p\left(x\right)=\exp\left(-\frac{1}{2}\left(x-\mu\right)^{\top}P^{-1}\left(x-\mu\right)\right)

    and measurement error are assumed to obey conditional Gaussian PDF

    \displaystyle p\left(d|x\right)=\exp\left(-\frac{1}{2}\left(d-Hx\right)^{\top}R^{-1}\left(d-Hx\right)\right)

    then the posterior PDF is given by

    \displaystyle p\left(x|d\right)\propto p\left(x\right)p\left(d|x\right)

    exponent for the posterior is

    \displaystyle -2E=x^{\top}\left(P^{-1}+H^{\top}R^{-1}H\right)x+x^{\top}\left(P^{-1}\mu+H^{\top}R^{-1}d\right)-\left(\mu^{\top}P^{-1}+d^{\top}R^{-1}H\right)x+\mu^{\top}P{}^{-1}\mu+d^{\top}R^{-1}d

    comparing quadratic terms in (posterior) PDF and following equivalent

    \displaystyle -2E=x^{\top}\hat{P}^{-1}x-x^{\top}\hat{P}^{-1}\hat{\mu}-\hat{\mu}^{\top}\hat{P}^{-1}x+\hat{\mu}^{\top}\hat{P}^{-1}\hat{\mu}

    posterior covariance is obtained as

    \displaystyle \hat{P}^{-1}=\left(P^{-1}+H^{\top}R^{-1}H\right)

    and comparing linear terms posterior mean is obtained as

    \displaystyle \hat{P}^{-1}\hat{\mu}=\left(P^{-1}\mu+H^{\top}R^{-1}d\right)

    or \hat{\mu}=\hat{P}\left(P^{-1}\mu+H^{\top}R^{-1}d\right).

    After some manipulation posterior covariance can be expresses as

    \hat{P}=\left(I-KH\right)P

    where K=PH^{\top}\left(R+HPH^{\top}\right)^{-1}

    and posterior mean can be expressed as

    \hat{\mu}=\mu+K\left(d-H\mu\right)

    Animation below illustrates evolution of updated state-error covariance matrix (light green) \hat{P}^\prime as the one of the principal standard deviations \sigma_2(\hat{P}) of prior state error covariance matrix (dark green) \hat{P} is varied. As expected, resulting (updated) covariance matrix approaches to state-estimate covariance matrix due to three measurements \hat{P} (orange ellipse). Also note that initially, when \sigma_2(\hat{P})=0, resulting state update resides along the straight line of initial state error covariance matrix (dark green) \hat{P} and along the dominant direction of state-estimate covariance matrix due to three measurements \hat{P} (orange ellipse).

    State update with three measurements as the prior state-error covariance matrix is varied

    Case of Measurements alone

    Measurements alone lead to an estimate of state with covariance

    P_{m}^{-1}=H^{\top}R^{-1}H

    and mean

    \hat{x}=P_{m}H^{\top}R^{-1}z

    Subscript “m” stands for “measurement”. Note the similarity to the least squares solution

    \hat{x}_{LS}=\left(H^{\top}H\right)^{-1}Hz and

    weighted LS solution

    \hat{x}_{WLS}=\left(H^{\top}WH\right)^{-1}H^{\top}Wz

    animation below illustrates estimated state \hat{x} from three measurement as the third (blue) measurement variance is increased. As expected mean (estimated) state goes away from blue line (measurement) to goes near other (green and red) lines (measurements). Also note that initially covariance ellipse is aligned with line of third measurement (blue) and as the standard deviation increases covariance ellipse becomes more circular with decrease in alignment/agreement with 3rd measurement.

    State estimate and covariance With 3 measurements with variance of a measurement (3rd~blue) varying

    Basic Examples of Kalman Filter

    State Error Covariance Evolution for Pendulum

    \displaystyle \mathbf{\dot{X}}=\mathbf{F}\left(\mathbf{X}\right)=\left[\begin{array}{c} \ddot{\theta}\\ u-\Omega^{2}\sin\theta\end{array}\right] = \left[\begin{array}{c} X_{2}\\ u-\Omega^{2}\sin X_{1}\end{array}\right]

    %% script for generating state trajectories of pendulum with covariance
    %{
    X0          :  initial state (nominal)
    nsamples    : number of samples in state point cloud (for covariance P0 and mean X0)
     x         :    cloud points (initial)
     P0         :   state-error covariance (initial)
    Tnom        :   number of steps nominal state trajectory
    Xnom        :   nominal trajectory states
    R1,R2       :  , measurement error covariance matrices
    S1,S2       :  , measurement only covariance estimates in state-space
    z1nom,z2nom :  , nominal trajectory measurements (nominal)
    v1,v2       :  , measurement noises for 1st and 2nd measurement models
    z1,z2       :  , measurements for 1st and 2nd models (z1nom/z2nom disturbed by v1/v2)
    Tcloud      :   time vectors for cloud trajectories
    Xcloud      :    state vectors for cloud trajectories
    %}
    
    global Omega
    
    alpha(0.5);
    
    Omega = 3;
    
    %% render state-derivative field plot
    mx1 = 12;
    mx2 = 12;
    [X1,X2] = meshgrid( linspace(-3,3,mx1), linspace(-3,3,mx2) );
    
    X1f = X1(:);
    X2f = X2(:);
    X12 = [X1f,X2f]';
    
    DX12 = pendulum(0,X12,0);
    
    phasePlanePlot = figure;
    hold on;
    quiver( X1f, X2f, DX12(1,:)', DX12(2,:)' );
    
    %% define initial covariance via its principal parameters
    %X0 = [-2;-1];
    X0 = [-2;-1]/3;
    sigma1b = 0.3/1;
    sigma2b = 0.2/1;
    thetab = 45*pi/180;
    P0 = generateCov1(sigma1b,sigma2b,thetab);
    
    [s1check, s2check, thetacheck] = analyseCov(P0);
    disp([s1check^2,s2check^2])
    
    %% generate state point cloud with Gaussian PDF of covariance P0
    nsamples = 54;
    randn('seed',10); %#ok
    x = generateMVG(P0,nsamples);
    
    % displace zero-mean points to initial state X0
    x(1,:) = x(1,:)+X0(1);
    x(2,:) = x(2,:)+X0(2);
    
    % check covariance estimation of point cloud
    P0est = cov(x);
    [sigma1bEst,sigma2bEst,thetabEst] = analyseCov(P0est);
    
    %% declare animated graphical entities (whose x,y,z data fields are modified to animate
    % function h = drawEllipse(c,r1,r2,theta,m)
    h = drawEllipse(X0, 3*sigma1b, 3*sigma2b, thetab, 32, []);
    hEst = drawEllipse(X0, 3*sigma1bEst, 3*sigma2bEst, thetabEst, 32, []);
    set(h,'Color','r','linewidth',3);
    
    measurementLine1 = line('color','c','linewidth',2,'linestyle','-.');
    measurementLine2 = line('color','m','linewidth',2,'linestyle','-.');
    measurementLine3 = line('color','y','linewidth',2,'linestyle','-.');
    
    measurementMean = plot(0,0,'^k','markerfacecolor','k','markersize',20);
    
    samplesplot = plot(x(1,:),x(2,:),'.','markersize',24); %
    meanplot = plot(X0(1),X0(2),'pr','markersize',25,'markerfacecolor','m');
    
    %% compute nominal trajectory "Xnom" for time range "Tnom" and plot time history
    Tcommon = linspace(0,2,4*128)';
    [Tnom,Xnom] = ode45(@pendulum, Tcommon, X0, odeset(), 0);
    
    timeHistoryPlot = figure;
    hold on;
    plot(Tnom,Xnom(:,1), Tnom, Xnom(:,2));
    
    % -------------------------------------------------
    % declare two measurement models to experiment with
    H1 = [[1,2];
        [3,-2]];
    
    H2 = [[1,2];
        [3,-2];
        [-1,1]];
    
    dt = mean(diff(Tnom));
    
    % -------------------------------------------------
    % define covariance for measurement error
    R1 = generateCov1(0.2,0.3,0*pi/180);
    R2 = diag([0.1,0.2,0.3])/4;
    
    S1 = inv(H1'*(R1\H1)); % inv(H1'*inv(R1)*H1);
    S2 = inv(H2'*(R2\H2));
    
    % -------------------------------------------------
    % compute measurements
    z1nom = H1*Xnom';
    z2nom = H2*Xnom';
    
    % -------------------------------------------------
    % create noise for measurements with covariance R1 and R2
    v1 = generateMVG( R1, length(Tcommon) );
    v2 = generateMVG( R2, length(Tcommon) );
    
    % -------------------------------------------------
    % add noise to measurements
    z1 = z1nom+v1;
    z2 = z2nom+v2;
    
    % -------------------------------------------------
    % plot nominal trajectory in state-space
    figure(phasePlanePlot);
    plot(Xnom(:,1),Xnom(:,2),'g','linewidth',4);
    
    % -------------------------------------------------
    % allocate space for cloud point trajectories (using repmat)
    Tcloud = repmat(Tcommon, [1,nsamples]);
    Xcloud = repmat(Xnom, [1,1,nsamples]);
    
    % -------------------------------------------------
    % simulate each cloud point and store in Xcloud along third dimension Xcloud(:,:,#cloudPoint)
    for i=1:size(x,2)
        % for all cloud points do
        [Tcloud(:,i),Xcloud(:,:,i)] = ode45(@pendulum, Tcommon, x(:,i), odeset(), 0);
        if mod(i,3)==0
            figure(phasePlanePlot);
            plot(Xcloud(1:3:end,1,i),Xcloud(1:3:end,2,i),'m:','linewidth',1);
        end;
        if 0
            figure(timeHistoryPlot);
            plot(Tcloud(:,i), Xcloud(:,1,i), Tcloud(:,i), Xcloud(:,2,i));
        end;
    end;
    
    % compute Jacobians at nominal trajectories
    [dxnom, dxxnom, dxunom] = pendulum(0, Xnom', 0);
    
    %% main animation loop (also information gathering for Asymptote animation)
    measurementModel = 1;
    
    % Asymptote rendering parameters pre-allocation
    allPi = zeros(length(Tcommon),3);
    allmi = zeros(length(Tcommon),2);
    allPiEst = zeros(length(Tcommon),3);
    alleta1 = zeros(2,length(Tcommon));
    alleta2 = zeros(2,length(Tcommon));
    
    figure(phasePlanePlot);
    for i=1:length(Tcommon)
        % take all cloud point coordinates ~ (x1i,x2i) at current time
        x1i = squeeze(Xcloud(i,1,:));
        x2i = squeeze(Xcloud(i,2,:));
        x12i = [x1i,x2i];
        % determine covariance matrix for cloud points at current time
        Pi = cov(x12i);
        mi = mean(x12i);
    
        allmi(i,:) = mi;
    
        % determine Jacobians for the current step
        Jx = dxxnom(:,:,i);
        Ju = dxunom(:,:,i);
        STM = expm(dt*Jx);
        PiEst = STM*P0*STM';
        P0 = PiEst;
    
        % analyse computed covariance -> principal components and rotation
        [s1i,s2i,thetai] = analyseCov(Pi);
        [s1iEst,s2iEst,thetaiEst] = analyseCov(PiEst);
    
        allPi(i,:) = [s1i,s2i,thetai];
        allPiEst(i,:) = [s1iEst,s2iEst,thetaiEst];
    
        % mean of measurements alone
        eta1 = S1*H1'*(R1\z1(:,i)); % S1*H1'*inv(R1)*z1(:,i);
        eta2 = S2*H2'*(R2\z2(:,i)); % S2*H2'*inv(R2)*z2(:,i);
    
        alleta1(:,i) = eta1;
        alleta2(:,i) = eta2;
    
        % update cloud points in drawing (for animation)
        set(samplesplot,'xdata',x1i);
        set(samplesplot,'ydata',x2i);
        set(samplesplot,'color',[0.4 0.5 0.7]);
    
        % update covariance matrix ellipsoid
        drawEllipse(Xnom(i,:), 3*s1i, 3*s2i, thetai, 32, h);
        drawEllipse(Xnom(i,:), 3*s1iEst, 3*s2iEst, thetaiEst, 32, hEst);
    
        % update estimated (mean) point
        set(meanplot, 'xdata', Xnom(i,1));
        set(meanplot, 'ydata', Xnom(i,2));
    
        switch measurementModel
            case 1
                set(measurementMean,'xdata', eta1(1));
                set(measurementMean,'ydata', eta1(2));
            case 2
                set(measurementMean,'xdata', eta2(1));
                set(measurementMean,'ydata', eta2(2));
        end;
    
        % draw computed mean
        plot(mi(1),mi(2),'.');
    
        %drawImplicitLine(drawImplicitLineNearThePoint(a,b,c,R,l,lh)
        switch measurementModel
            case 1
                drawImplicitLineNearThePoint(H1(1,1), H1(1,2), z1(1,i), Xnom(i,:)', 3, measurementLine1);
                drawImplicitLineNearThePoint(H1(2,1), H1(2,2), z1(2,i), Xnom(i,:)', 3, measurementLine2);
            case 2
                drawImplicitLineNearThePoint(H2(1,1), H2(1,2), z2(1,i), Xnom(i,:)', 3, measurementLine1);
                drawImplicitLineNearThePoint(H2(2,1), H2(2,2), z2(2,i), Xnom(i,:)', 3, measurementLine2);
                drawImplicitLineNearThePoint(H2(3,1), H2(3,2), z2(3,i), Xnom(i,:)', 3, measurementLine3);
        end;
    
        if 1
            axis([-4,4,-4,4]);
        else
            axis( [Xnom(i,1)-2 Xnom(i,1)+2 Xnom(i,2)-2 Xnom(i,2)+2] );
        end;
    
        drawnow;
        %input('>>');
    end;
    
    %% Asymptote code generation
    
    allCloudX1 = squeeze(Xcloud(:,1,:));
    allCloudX2 = squeeze(Xcloud(:,2,:));
    
    fid = fopen('asycode.asy','wt');
    
    fprintf(fid,'real allPi[][] = \n');
    fprintf(fid,'%s;\n',Matlab2AsymptoteMatrix(allPi));
    fprintf(fid,'/* ____________________________ */\n');
    
    fprintf(fid,'real allPiEst[][] = \n');
    fprintf(fid,'%s;\n',Matlab2AsymptoteMatrix(allPiEst));
    fprintf(fid,'/* ____________________________ */\n');
    
    fprintf(fid,'real allmi[][] = \n');
    fprintf(fid,'%s;\n',Matlab2AsymptoteMatrix(allmi));
    fprintf(fid,'/* ____________________________ */\n');
    
    fprintf(fid,'real alleta1[][] = \n');
    fprintf(fid,'%s;\n',Matlab2AsymptoteMatrix(alleta1 ));
    fprintf(fid,'real alleta2[][] = \n');
    fprintf(fid,'%s;\n',Matlab2AsymptoteMatrix(alleta2));
    fprintf(fid,'/* ____________________________ */\n');
    
    fprintf(fid,'real Xnom[][] = \n');
    fprintf(fid,'%s;\n',Matlab2AsymptoteMatrix(Xnom));
    fprintf(fid,'/* ____________________________ */\n');
    
    fprintf(fid,'real cloudX1[][] = \n');
    fprintf(fid,'%s;\n',Matlab2AsymptoteMatrix(allCloudX1));
    fprintf(fid,'real cloudX2[][] = \n');
    fprintf(fid,'%s;\n',Matlab2AsymptoteMatrix(allCloudX2));
    fprintf(fid,'/* ____________________________ */\n');
    
    fclose(fid);
    

    BayesianApproachToKF_save2

    Mechanism Analysis and Animation

    January 25, 2010

    This post illustrates a mechanism analysis technique, known as “Modulmetode”[1]. In this method aim is to construct unknown points (joint positions) starting from the known ones (driven linkages/joints, fixed joints) and applying a formula available for the motion of sub-assemblies of the mechanism. [1] gives a schematic list of some kinematic pairs (or motion types) which are used as basic entities in “Modulmetode”. In this post only one of the schemas (“Zweischlag”) is illustrated.

    “Zweischlag” modul

    Linkages in this kinematic pair are connected by a revolute joint. Also position of distinct fixed points for each of the linkages together with distance to common (revolute joint) are defined. Sought variable is the common revolute joint position. Figure below depicts this problem and then follows the solution.

    third point
    Position of unknown point is given relative to “base” segment that connects the known positions.

    \displaystyle \alpha = \frac{x_{1}}{d} =\frac{d^{2}+d_{1}^{2}-d_{2}^{2}}{2d^{2}}=\frac{d^{2}+\left(  d_{1}-d_{2}\right)  \left(  d_{1}+d_{2}\right)  }{2d^{2}}

    \displaystyle \beta = \frac{h}{d}  =\frac{1}{2d^{2}}\sqrt{\left(  \left(  d_{2}+d_{1}\right)  ^{2} -d^{2}\right)  \left(  d^{2}-\left(  d_{2}-d_{1}\right)  ^{2}\right)  }

    For auxiliary variables defined by

    \displaystyle \frac{d_{1}-d_{2}}{d}  =\delta

    \displaystyle \frac{d_{1}+d_{2}}{d}  =\sigma

    local normalized coordinates can be represented by

    \displaystyle \alpha =\frac{1+\sigma\delta}{2}

    \displaystyle \beta =\frac{1}{2}\sqrt{\left(  \sigma^{2}-1\right)  \left(  1-\delta^{2}\right)  }

    Then for given positions of driven joints of two connected linkages z_i,z_o, common joint locations is given by

    \eta\left(  z_{1},z_{2},d_{1},d_{2}\right)  =z_{1}+\left(  \alpha\pm i\beta\right)  \left(  z_{2}-z_{1}\right)

    For the fourbar mechanism analysis is simply

    \displaystyle P_{i}    =G_{i}+ae^{i\theta}

    \displaystyle G_{o}   =G_{i}+d

    \displaystyle P_{o}   =\eta\left(  P_{i},G_{o},b,c\right)

     

     

    Asymptote source code to analyze and animate fourbar is given below

    /* FILE fourbarMP1.asy */
    import animate;
    import patterns;     // for coupler body
    import roundedpath;  // for coupler body
    import geometry;		// for distance()
    
    bool multipage = false;
    
    settings.tex="pdflatex";
    settings.outformat="pdf";
    
    add("crosshatch",crosshatch(3mm));
    
    if(multipage){
    	settings.keep=true;
    } else {
    	settings.keep=true;
    }
    
    animation anim;
    
    pair rates(pair z1, pair z1dot, pair z2, pair z2dot){
    	/* returns the distance and rate of distance between points z1 and z2 */
    	pair delta = z1-z2;
    	pair ddelta = z1dot-z2dot;
    	real d = abs(z1-z2);
    	real ddot = dot(delta,ddelta)/d;
    	return (d,ddot);
    };
    
    pair[] tpaux(real d1, real d2, real d, real ddot, real sign){
    	/*
    	In 1st argument returns normalized (complex) coordinate
    	of point at a distance d1 and d2 from two points that
    	are apart by distance d. In 2nd argument returns the rate
    	of first argument. Upper/lower configuration of solution
    	is selected by input argument sign=-1,+1
    	*/
    
    	real sigma = (d1+d2)/d;
    	real delta = (d1-d2)/d;
    
    	real alpha = (1+delta*sigma)/2;
    	real beta = (1-delta**2)*(sigma**2-1);
    
    	if(beta<0){
    		beta = 0.0;
    	} else {
    		beta = 0.5*sqrt( beta );
    	}
    
    	real alphadot = -sigma*delta*(ddot/d);
    	real betadot = (1./(4.0*beta))*(sigma^2+delta^2-2*sigma^2*delta^2)*(ddot/d);
    
    	pair[] result = {(alpha,beta*sign),(alphadot,betadot*sign)};
    
    	return result;
    }
    
    pair[] FourbarJointsAndCouplerPoint(real l1, real l2, real l3, real l4, real theta, real thetadot, pair CP, real sign){
    	/*
    	In first 4 arguments returns the joint positions of fourbar.
    	In 5th argument returns the coupler point, described
    	coupler relative (complex) coordinate CP.
    	In the 6th argument returns the projection of POUT to line G2-PIN.
    	*/
    
    	pair G1 = (0,0);
    	pair PIN = (l1*Cos(theta),l1*Sin(theta));
    	pair dPIN = (-l1*Sin(theta),l1*Cos(theta))*thetadot;
    	pair G2 = (l4,0);
    
    	pair ddd = rates(PIN,dPIN,G2,(0,0));  // without rate this would be simply abs(PIN-G2);
    	real d1 = l2;
    	real d2 = l3;
    	pair[] res = tpaux(d1,d2,ddd.x,ddd.y,sign);
    	pair ab = res[0];
    	pair abdot = res[1];
    
    	pair POUT = PIN + ab * (G2-PIN);
    
    	pair[] result = {G1 , PIN , POUT , G2 , PIN+CP*(POUT-PIN), PIN + ab.x*(G2-PIN) };
    
    	return  result;
    }
    
    /* Define link lengths */
    real LL1,LL2,LL3,LL4;
    LL1 = 0.75;
    LL2 = 2;
    LL3 = 4;
    LL4 = 3;
    
    /* Set number of frames */
    int nof=25;
    
    /* Define a coupler point (normalized and relative to coupler) */
    pair CP = (1.,0.55);
    
    /* Define the solution configuration */
    real sign = +1;
    
    path linkages;
    real theta;  /* input link rotation angle */
    guide couplerPoints;
    path couplerBody;
    
    picture pic;
    size(pic,10cm);
    
    for(int i=0; i < nof; ++i) {
      write("rendering frame: ", i);
    
      theta = i*(360./(nof-1));
    
      /* Determine fourbar joint positions and coupler point */
      pair g[]=FourbarJointsAndCouplerPoint(LL1,LL2,LL3,LL4,theta,1.,CP,sign);
    
      couplerPoints = couplerPoints -- g[4];
    
      /* set linkages path */
      linkages = g[0] -- g[1] -- g[2] -- g[3];
    
      /* draw fourbar and coupler point and frame */
      draw(pic,linkages,red+5bp);			// draw linkages as lines
      dot(pic,linkages,heavyred+10bp);	// draw linkage joints as points
      label(pic, "$G_0$", g[0], 2W, deepcyan+fontsize(8));
      label(pic, "$P_I$", g[1], 2N, deepcyan+fontsize(8));
      label(pic, "$P_O=P_I \pm \left(\alpha+i\beta\right)(G_1-P_I)$", g[2], E+N, deepcyan+fontsize(8));
      label(pic, "$G_1$", g[3], 2E, deepcyan+fontsize(8));
      dot(pic,"$CP$",g[4],2N,magenta+8bp);	// draw coupler point as point
      draw(pic, couplerPoints, cmyk(green) );
    
      if(false){
      	draw(pic, Label("$\alpha d$", position=MidPoint), g[1] -- g[5], dashed+gray+black);
      	draw(pic, Label("$\left(1-\alpha\right) d$", position=MidPoint), g[5] -- g[3], dashed+gray+black);
      	draw(pic, Label("$\beta d$"), g[5] -- g[2], dashed+gray+black);
      } else {
      	draw(pic, g[1] -- g[5], dashed+gray+black);
      	draw(pic, g[5] -- g[3], dashed+gray+black);
      	draw(pic, g[5] -- g[2], dashed+gray+black);
    
      	distance(pic, "$\alpha d$", g[1], g[5], -4mm, gray);  // +fontsize(8)
      	distance(pic, "$\left(1-\alpha\right) d$", g[5], g[3], 6mm, gray); //+fontsize(8)
      	distance(pic, "$\beta d$", g[5], g[2], 4mm, gray); //+fontsize(8)
      }
    
      couplerBody = roundedpath(g[1] -- g[2] -- g[4] -- cycle, 0.2);
      filldraw( pic, couplerBody, royalblue+white+opacity(0.6), royalblue+1bp );  // pattern("crosshatch")
    
      anim.add(pic);
    
      erase(pic);
    }
    if(multipage){
    	anim.pdf("",multipage=true);
    } else {
    	label(anim.pdf(BBox(1mm,nullpen),delay=500,"controls,loop"));
    }
    

    when executed this file generates

  • fourbarMP1.pdf (since multipage=false)
  • _fourbarMP1.pdf (since settings.keep=true)
  • other intermediate files _fourbarMP1_0.pdf to _fourbarMP1_24.pdf and _fourbarMP1_0.eps to _fourbarMP1_24.eps
  • To create a gif animation from intermediate files, first convert intermediate PDF files to some bitmap format, say PNG, by ImageMagick

    rem convert each PDF to PNG
    convert -filter Mitchell -resize 400x400 -density 300x300 _fourbarMP1_%d.pdf[0-24] FB_%d.png
    
    rem convert/combine generated PNG files to animated GIF
    convert FB_%d.png[0-24] -loop 0 gifanim.gif
    

    In this case, strangely, intermediate PDF files do not contain labels, but the final PDF animation files do. For this reason, to get correct output, PNG files are created from multipage PDF file (ImageMagick can automatically generate a bitmap file for each page of multipage pdf). Correction to above command line is then

    rem generate PNG for each page of multipage PDF
    convert -filter Mitchell -resize 400x400 -density 300x300 _fourbarMP1.pdf _fourbarMP1_%02d.png
    
    rem convert/combine generated PNG to animated GIF
    convert _fourbarMP1_%02d.png[0-24] -loop 0 gifanim.gif
    

    Resulting animated GIF file is given below.

    Animation of fourbar and coupler point

     

    Note that Asymptote already seems to have a mechanism to pass arbitrary command line switches to ImageMagick by “convertOptions”. But using this option would not be as convenient as the framework given above if the rendering time is long. Manual conversion of pdf/eps files allows quick experimentation of conversion factors (quality, size, etc).

    PDF versions of above animation are _fourbarMP1 (Multi-page) and fourbarMP1 (Single-page)

    Differential Kinematics in “Modulmetode”

    To compute the rate of output point, rates of auxiliary variables, given below, are required

    \displaystyle \dot{\delta}   =-\frac{d_{1}-d_{2}}{d^{2}}\dot{d}=-\delta\frac{\dot{d}}{d}

    \displaystyle \dot{\sigma}   =-\frac{d_{1}+d_{2}}{d^{2}}\dot{d}=-\sigma\frac{\dot{d}}{d}

    \displaystyle \dot{\alpha} =-\sigma\delta\frac{\dot{d}}{d}

    \displaystyle \dot{\beta} =-\frac{1}{4\beta}\frac{\dot{d}}{d}\left(  \sigma^{2}-2\sigma^{2}\delta^{2}+\delta^{2}\right)

     

    Polode in “Modulmetode”

    If a point on coupler link is described by z=z_{i}+\rho\left(  z_{o}-z_{i}\right) then \dot{z}=\dot{z}_{i}+\rho\left(  \dot{z}_{o}-\dot{z}_{i}\right) gives the polode location (coupler frame relative) as

    \displaystyle \rho_{0}=-\frac{\dot{z}_{i}}{\dot{z}_{o}-\dot{z}_{i}}

    and in fixed frame

    \displaystyle z_{0} =z_{i}+\frac{\dot{z}_{i}}{\dot{z}_{o}-\dot{z}_{i}}\left(  z_{i}-z_{o}\right)

    In general dyadic/matrix formulation[2], for a motion given by

    \mathbf{P}=\mathbf{Ap}+\mathbf{d}

    differential kinematics

    \mathbf{\dot{P}} =\mathbf{\dot{A}p+\dot{d}}

    yields the polode position (instantaneous rotation center in fixed space)

    \mathbf{P}^{\ast}=\mathbf{d}-\mathbf{\Omega}^{-1}\mathbf{\dot{d}}

    and the centrode position (instantaneous rotation center in moving space)

    \mathbf{p}^{\ast}  =-\left(  \Omega\mathbf{A}\right)  ^{-1}\mathbf{\dot{d}}

    animation below shows centrode rolling over (without slipping) polhode for a fourbar with link dimensions a = 2.25, b = 2.8, c = 4, d= 1

    Centrode rolling over Polhode

    Note that in literature, instead of terms polhode and centrode, fixed polhode/fixed centrode and moving polhode/moving centrode may be encountered. The meaning of these terms in this post conforms to that in [2].

    Source Code To Render Polode Animation

    Following code defines routine to draw connecting rod schematically

    import graph;  // for Arc()
    
    size(10cm);
    
    real r1,r2,R1,R2,rho1,rho2,d;
    
    r1 =  0.2;
    r2 =  0.4;
    R1 =  0.3;
    R2 =  0.5;
    rho1 = 0.3;
    rho2 = 0.3;
    d = 1;
    
    real sqr(real x){
    	return x*x;
    };
    
    void drawCenterSign(pair center, real angle, real rad, real ext, pen pen1, pen pen2, pen pen3){
    	transform t = shift(center)*rotate(angle);
    
    	fill( t*( (0,0) -- arc( (0,0), rad,   0,  90) -- cycle ) , pen1);
    	fill( t*( (0,0) -- arc( (0,0), rad,  90, 180) -- cycle ), pen2);
    	fill( t*( (0,0) -- arc( (0,0), rad, 180, 270) -- cycle ), pen1);
    	fill( t*( (0,0) -- arc( (0,0), rad, 270, 360) -- cycle ), pen2);
    
    	real L = rad+ext;
    	draw( t*( (0,L) -- (0,-L)), pen3);
    	draw( t*( (L,0) -- (-L,0)), pen3);
    
    }
    
    path[] ConnectingRod(pair origin, real angle, real r1, real r2, real R1, real R2, real rho1, real rho2, real d, real hole1, real hole2){
    	transform t = shift(origin)*rotate(angle);
    	real a = asin((r2-r1)/d);
    	pair p1 = r1*exp(I*(pi/2+a));
    	pair p2 = d+r2*exp(I*(pi/2+a));
    
    	pair q1 = p1 + sqrt(sqr(R1)-sqr(r1))*exp(I*a);
    	pair q2 = p2 - sqrt(sqr(R2)-sqr(r2))*exp(I*a);
    
    	real p1s1 = sqrt( sqr(R1+rho1)-sqr(rho1+r1) );
    	real p2s2 = sqrt( sqr(R2+rho2)-sqr(rho2+r2) );
    
    	real b1 = acos((rho1+r1)/(rho1+R1));
    	real b2 = acos((rho2+r2)/(rho2+R2));
    
    	pair w1 = (rho1+R1)*exp(I*(pi/2+a-b1));
    	pair w2 = (d,0) + (rho2+R2)*exp(I*(pi/2+a+b2));
    
    	pair s1 = p1 + p1s1*exp(I*a);
    	pair s2 = p2 - p2s2*exp(I*a);
    
    	path[] segments = new path[8+2]; // +2 for holes
    
    	if(false){
    		segments[0] = arc( (0,0), R1, degrees(2*pi-(pi/2+a-b1)), degrees(pi/2+a-b1));
    		segments[1] = arc( w1   , rho1, degrees(3*pi/2+a-b1) , degrees(3*pi/2+a));
    		segments[2] = s1 -- s2;
    		segments[3] = arc( w2,    rho2, degrees(3*pi/2+a) , degrees(3*pi/2+a+b2));
    		segments[4] = arc( (d,0), R2, degrees(pi/2+a+b2), degrees(2*pi-(pi/2+a+b2)), CW);
    
    		segments[5] = arc( (w2.x,-w2.y),    rho2, degrees(2*pi-(3*pi/2+a)) , degrees(2*pi-(3*pi/2+a+b2)));
    		segments[6] = (s1.x,-s1.y) -- (s2.x,-s2.y);
    		segments[7] = arc( (w1.x, -w1.y)   , rho1, degrees(2*pi-(3*pi/2+a-b1)) , degrees(2*pi-(3*pi/2+a)));
    
    		if(hole1>0){
    			segments[8] = arc( (0,0), hole1, 0, 360);
    		} else {
    			segments[8] = nullpath;
    		}
    		if(hole2>0){
    			segments[9] = arc( (d,0), hole2, 0, 360);
    		} else {
    			segments[9] = nullpath;
    		}
    
    	} else {
    		segments[0] = Arc( (0,0), R1, degrees(2*pi-(pi/2+a-b1)), degrees(pi/2+a-b1));
    		segments[1] = Arc( w1   , rho1, degrees(3*pi/2+a-b1) , degrees(3*pi/2+a));
    		segments[2] = s1 -- s2;
    		segments[3] = Arc( w2,    rho2, degrees(3*pi/2+a) , degrees(3*pi/2+a+b2));
    		segments[4] = Arc( (d,0), R2, degrees(pi/2+a+b2), degrees(2*pi-(pi/2+a+b2)), CW);
    
    		segments[5] = Arc( (w2.x,-w2.y),    rho2, degrees(2*pi-(3*pi/2+a)) , degrees(2*pi-(3*pi/2+a+b2)));
    		segments[6] = (s1.x,-s1.y) -- (s2.x,-s2.y);
    		segments[7] = Arc( (w1.x, -w1.y)   , rho1, degrees(2*pi-(3*pi/2+a-b1)) , degrees(2*pi-(3*pi/2+a)));
    
    		if(hole1>0){
    			segments[8] = Arc( (0,0), r1, 0, 360);
    		} else {
    			segments[8] = nullpath;
    		}
    		if(hole2>0){
    			segments[9] = Arc( (d,0), r2, 0, 360);
    		} else {
    			segments[9] = nullpath;
    		}
    
    	}
    
    	for(int i=0;i<segments.length;++i){
    		segments[i] = t*segments[i];
    	}
    
    	return segments;
    };
    
    void drawConnectingRod(pair origin, real angle, real r1, real r2, real R1, real R2, real rho1, real rho2, real d,bool singlePath, real hole1, real hole2, pen[] p){
    	path[] conrodshape = ConnectingRod(origin, angle, r1, r2, R1, R2, rho1, rho2, d, hole1, hole2);
    	if(singlePath){
    		draw( conrodshape,  p[0] );
    	} else {
    
    		draw( conrodshape[0], p[0]);
    		draw( conrodshape[1], p[1]);
    		draw( conrodshape[2], p[2]);
    		draw( conrodshape[3], p[3]);
    		draw( conrodshape[4], p[4]);
    		draw( conrodshape[5], p[5]);
    		draw( conrodshape[6], p[6]);
    		draw( conrodshape[7], p[7]);
    		draw( conrodshape[8], p[8]);
    		draw( conrodshape[9], p[9]);
    	}
    }
    
    bool singlePath = false;
    pen basicPen = linewidth(3bp);
    pen[] pens={ basicPen+red,
                 basicPen+blue,
                 basicPen+magenta,
                 basicPen+blue,
                 basicPen+red,
    
                 basicPen+green,
                 basicPen+yellow,
                 basicPen+green,
                 basicPen+royalblue,
                 basicPen+royalblue};
    
    drawConnectingRod(origin=(0,0), angle=0, r1, r2, R1, R2, rho1, rho2, d, singlePath, r1, r2, pens);
    drawCenterSign( center=(0,0), angle=0, rad=0.2, ext=0.05, black, yellow, black+1pt);
    
    drawConnectingRod(origin=(2,0), angle=20, r1, r2, R1, R2, rho1, rho2, d, singlePath, r1, -r2, pens);
    drawCenterSign( center=(2,0), angle=20, rad=0.2, ext=0.05, black, yellow, black+1pt);
    
    drawConnectingRod(origin=(2,2), angle=60, r1, r2, R1, R2, rho1, rho2, d, singlePath, -r1, r2, pens);
    drawCenterSign( center=(2,2), angle=60, rad=0.2, ext=0.05, black, yellow, black+1pt);
    

    This Asymptote code results in  this pdf file and may be converted to png file (by ImageMagick as above, for example)

    sample test call to connecting rod drawing routine

    schematic dimensions are illustrated with following Asymptote code

    import settings;
    import graph;  // for Arc()
    import geometry; // for markangle, dimension
    
    tex="latex";
    outformat = "pdf";
    keep=true;
    keepaux=true;
    
    size(15cm,0);
    unitsize(1cm);
    
    real r1,r2,R1,R2,rho1,rho2,d,sc;
    
    r1 =  2;
    r2 =  4;
    R1 =  3;
    R2 =  5;
    rho1 = 3;
    rho2 = 3;
    d = 12;
    
    real sqr(real x){
    	return x*x;
    };
    
    pair origin = (0,0);
    real angle = 0;
    bool singlePath = false;
    pen dimPen = 0.25pt+dashed;
    pen basicPen = linewidth(4bp);
    pen[] pens={ basicPen+red,
                 basicPen+blue,
                 basicPen+heavygreen,
                 basicPen+blue,
                 basicPen+red,
    
                 basicPen+blue,
                 basicPen+heavygreen,
                 basicPen+blue };
    
    	transform t = shift(origin)*rotate(angle);
    	real a = asin((r2-r1)/d);
    	pair p1 = r1*exp(I*(pi/2+a));
    	pair p2 = d+r2*exp(I*(pi/2+a));
    
    	pair q1 = p1 + sqrt(sqr(R1)-sqr(r1))*exp(I*a);
    	pair q2 = p2 - sqrt(sqr(R2)-sqr(r2))*exp(I*a);
    
    	real p1s1 = sqrt( sqr(R1+rho1)-sqr(rho1+r1) );
    	real p2s2 = sqrt( sqr(R2+rho2)-sqr(rho2+r2) );
    
    	real b1 = acos((rho1+r1)/(rho1+R1));
    	real b2 = acos((rho2+r2)/(rho2+R2));
    
    	pair w1 = (rho1+R1)*exp(I*(pi/2+a-b1));
    	pair w2 = (d,0) + (rho2+R2)*exp(I*(pi/2+a+b2));
    
    	pair s1 = p1 + p1s1*exp(I*a);
    	pair s2 = p2 - p2s2*exp(I*a);
    
    	pair t1 = (R1*w1+rho1*(0,0))/(R1+rho1);
    	pair t2 = (R2*w2+rho2*(d,0))/(R2+rho2);
    
    	/*
    	arc(0,0,r1,r1,0,2*pi,'-.k',1);
    	arc(d,0,r2,r2,0,2*pi,'-.k',1);
    
    	arc(0,0,R1,R1,0,2*pi,'-.k',1);
    	arc(d,0,R2,R2,0,2*pi,'-.k',1);
    
    	arc(w1x,w1y,rho1,rho1,0,2*pi,'-.k',1);
    	arc(w2x,w2y,rho2,rho2,0,2*pi,'-.k',1);
    	*/
    
    	path[] segments = new path[8];
    
    	if(false){
    		segments[0] = arc( (0,0), R1, degrees(2*pi-(pi/2+a-b1)), degrees(pi/2+a-b1));
    		segments[1] = arc( w1   , rho1, degrees(3*pi/2+a-b1) , degrees(3*pi/2+a));
    		segments[2] = s1 -- s2;
    		segments[3] = arc( w2,    rho2, degrees(3*pi/2+a) , degrees(3*pi/2+a+b2));
    		segments[4] = arc( (d,0), R2, degrees(pi/2+a+b2), degrees(2*pi-(pi/2+a+b2)), CW);
    
    		segments[5] = arc( (w2.x,-w2.y),    rho2, degrees(2*pi-(3*pi/2+a)) , degrees(2*pi-(3*pi/2+a+b2)));
    		segments[6] = (s1.x,-s1.y) -- (s2.x,-s2.y);
    		segments[7] = arc( (w1.x, -w1.y)   , rho1, degrees(2*pi-(3*pi/2+a-b1)) , degrees(2*pi-(3*pi/2+a)));
    	} else {
    		segments[0] = Arc( (0,0), R1, degrees(2*pi-(pi/2+a-b1)), degrees(pi/2+a-b1));
    		segments[1] = Arc( w1   , rho1, degrees(3*pi/2+a-b1) , degrees(3*pi/2+a));
    		segments[2] = s1 -- s2;
    		segments[3] = Arc( w2,    rho2, degrees(3*pi/2+a) , degrees(3*pi/2+a+b2));
    		segments[4] = Arc( (d,0), R2, degrees(pi/2+a+b2), degrees(2*pi-(pi/2+a+b2)), CW);
    
    		segments[5] = Arc( (w2.x,-w2.y),    rho2, degrees(2*pi-(3*pi/2+a)) , degrees(2*pi-(3*pi/2+a+b2)));
    		segments[6] = (s1.x,-s1.y) -- (s2.x,-s2.y);
    		segments[7] = Arc( (w1.x, -w1.y)   , rho1, degrees(2*pi-(3*pi/2+a-b1)) , degrees(2*pi-(3*pi/2+a)));
    
    	}
    
    	draw( arc( (0,0), R1, 0, 360), dimPen+red);
    	draw( arc( (0,0), r1, 0, 360), dimPen+royalblue);
    	draw( arc( (d,0), R2, 0, 360), dimPen+red);
    	draw( arc( (d,0), r2, 0, 360), dimPen+royalblue);
    	draw( arc( w1, rho1, 0, 360), dimPen+blue);
    	draw( arc( w2, rho2, 0, 360), dimPen+blue);
    
    	dot( "$P_1$", p1, NW);
    	dot( "$P_2$", p2, NW);
    	dot( "$W_1$", w1, N);
    	dot( "$W_2$", w2, N);
    	dot( "$S_1$", s1, SE);
    	dot( "$S_2$", s2, SE);
    	dot( "$T_1$", t1, 1.5N);
    	dot( "$T_2$", t2, 1.5N);
    	dot( "$(0,0)$", (0,0), S);
    	dot( "$(d,0)$", (d,0), S);
    
    	coordsys R=cartesiansystem(O=5*exp(I*a), i=1.2*exp(I*a), j=1.2*I*exp(I*a));
    	show("$O$","$e^{i\alpha}$", "$i e^{i\alpha}$", R, xpen=invisible);
    
    	pair aux = d+(r2-r1)*I*exp(I*a);
    	draw( (0,0) -- aux, dimPen+royalblue );
    	markangle( "$\alpha$", (d,0),(0,0),aux,radius=65mm, ArcArrow);
    
    	if(true){
    		draw( Label("$\frac{\pi}{2}+\alpha$",EndPoint), w1 -- w1+1.5*exp(I*(a+pi/2)), gray(0.7)+fontsize(6), EndArrow());
    		draw( Label("$\frac{\pi}{2}+\alpha-\beta_1$",EndPoint), w1 -- w1+1.5*exp(I*(a+pi/2-b1)), gray(0.7)+fontsize(6), EndArrow());
    
    		draw( Label("$\frac{\pi}{2}+\alpha$",EndPoint), w2 -- w2+1.5*exp(I*(a+pi/2)), gray(0.7)+fontsize(6), EndArrow());
    		draw( Label("$\frac{\pi}{2}+\alpha+\beta_2$",EndPoint), w2 -- w2+1.5*exp(I*(a+pi/2+b2)), gray(0.7)+fontsize(6), EndArrow());
    	}
    
    	draw( (-R1,0) -- (d+R2,0), dimPen+royalblue);
    	distance( "$d$", (0,0), (d,0), 9mm, joinpen=dotted);
    
    	draw( t1--w1--s1, dimPen);
    	draw( t2--w2--s2, dimPen);
    	markangle( "$\beta_1$", t1,w1,s1, radius=10mm, ArcArrow);
    	markangle( "$\beta_2$", s2,w2,t2, radius=10mm, ArcArrow);
    
    	draw( (0,R1)--(0,0)--(R1*I*exp(I*a)), dimPen);
    	draw( (d,R2)--(d,0)--(d+R1*I*exp(I*a)), dimPen);
    	markangle( "$\alpha$", (0,R1),(0,0),p1, radius=18mm, ArcArrow);
    	markangle( "$\alpha$", (d,R2),(d,0),p2, radius=20mm, ArcArrow);
    
    	draw( w1 -- s1, dimPen);
    	distance( "$\rho_1$", w1, s1,-3mm,joinpen=dotted);
    	draw( w2 -- s2, dimPen);
    	distance( "$\rho_2$", w2, s2,+3mm,joinpen=dotted);
    
    	draw( s1 -- (0,0)+p1s1*exp(I*a), dimPen);
    	draw( s2 -- (d,0)-p2s2*exp(I*a) -- (d,0), dimPen);
    	draw( t1 -- (0,0), dimPen);
    	draw( t2 -- (d,0), dimPen);
    
    	draw( (0,0) -- p1, dimPen);
    	distance( "$r_1$", (0,0) , p1,-10mm,joinpen=dotted);
    	distance( "$R_1$", (0,0) , R1*I*exp(I*a),-30mm,joinpen=dotted);
    
    	distance( "$r_2$", (d,0) , p2,10mm,joinpen=dotted);
    	distance( "$R_2$", (d,0) , d+R2*I*exp(I*a),40mm,joinpen=dotted);
    
    	draw( p1 -- p2, dimPen+royalblue);
    
    	for(int i=0;i<segments.length;++i){
    		segments[i] = t*segments[i];
    	}
    
    	if(singlePath){
    		draw( segments, pens[0] );
    	} else {
    		draw( segments[0], pens[0]);
    		draw( segments[1], pens[1]);
    		draw( segments[2], pens[2]);
    		draw( segments[3], pens[3]);
    		draw( segments[4], pens[4]);
    		draw( segments[5], pens[5]);
    		draw( segments[6], pens[6]);
    		draw( segments[7], pens[7]);
    	}
    
    	draw( scale(0.85)* Label( minipage("\begin{eqnarray*}\sin \alpha  &=&\frac{r_{2}-r_{1}}{d} \\
    P_{1} &=&r_{1}e^{i\left( \frac{\pi }{2}+\alpha \right) }=ir_{1}e^{i\alpha }
    \\
    P_{2} &=&d+r_{2}e^{i\left( \frac{\pi }{2}+\alpha \right)
    }=d+ir_{2}e^{i\alpha } \\
    \cos \left( \beta _{1}\right)  &=&\frac{\rho _{1}+r_{1}}{\rho _{1}+R_{1}} \\
    \cos \left( \beta _{2}\right)  &=&\frac{\rho _{2}+r_{2}}{\rho _{2}+R_{2}} \\
    Q_{1} &=&P_{1}+\sqrt{R_{1}^{2}-r_{1}^{2}}e^{i\alpha } \\
    Q_{2} &=&P_{2}-\sqrt{R_{2}^{2}-r_{2}^{2}}e^{i\alpha } \\
    \overline{P_{1}S_{1}} &=&\sqrt{\left( R_{1}+\rho _{1}\right) ^{2}-\left(
    \rho _{1}+r_{1}\right) ^{2}} \\
    \overline{P_{2}S_{2}} &=&\sqrt{\left( R_{2}+\rho _{2}\right) ^{2}-\left(
    \rho _{2}+r_{2}\right) ^{2}} \\
    S_{1} &=&P_{1}+\overline{P_{1}S_{1}}e^{i\alpha } \\
    S_{2} &=&P_{2}+\overline{P_{2}S_{2}}e^{i\alpha } \\
    W_{1} &=&\left( R_{1}+\rho _{1}\right) e^{i\left( \frac{\pi }{2}+\alpha
    -\beta _{1}\right) } \\
    W_{2} &=&d+\left( R_{2}+\rho _{2}\right) e^{i\left( \frac{\pi }{2}+\alpha
    +\beta _{2}\right) } \\
    T_{1} &=&R_{1}e^{i\left( \frac{\pi }{2}+\alpha -\beta _{1}\right) } \\
    T_{2} &=&d+R_{2}e^{i\left( \frac{\pi }{2}+\alpha +\beta _{2}\right) }
    \end{eqnarray*}", 9.5cm), (6,-6), 2S, heavyblue), box, paleblue+2bp, xmargin=4mm,ymargin=2mm);
    

    which results in PDF file and converted to png

    connecting rod dimensions

     

    Implementing in Processing language

    Framework presented may be implemented in Processing language. A screenshot of Procesing IDE is given below.

    Processing IDE

    Sample implementation capable of producing on screen and pdf outputs is given below. Code can be executed on internet site openprocessing.org via copy paste to its editor. A saved version is stored in this location

    /* FILE: animateFourvar2.pde */
    import processing.pdf.*;
    
    /* global variables */
    // dimension of fourbar ... link lengths
    float a,b,c,d;
    // initial angle and angle steps of input crank angle
    float theta,dtheta;
    float[] IP,OP1,GP1,GP2,OP2;
    boolean pdfOutput = false;
    PGraphicsPDF pdf;
    
    void circularArc(float cx, float cy, float rad,float sa,float ea,float adelta){
      /*
      draw a circular arc centered at (cx,cy), having radius rad, angle sweep (sa,ea), and angle step adelta
      Supposed to be called within BeginShape .. EndShape
      */
      float newgx, newgy;
    
      int adiv = round( rad*abs(ea-sa)/adelta );
      adelta = (ea-sa)/adiv;
    
      float gx = rad*cos(sa);
      float gy = rad*sin(sa);
      float rx = cos(adelta);
      float ry = sin(adelta);
      float a = sa;
      for(int i=0;i<adiv;i++){
        if(false){
          // g <- (gx+I*gy)*(rx+I*ry) = (gx*rx-gy*ry)+I*(gx*ry+gy*rx)
          newgx = (gx*rx-gy*ry);
          newgy = (gx*ry+gy*rx);
          gx = newgx;
          gy = newgy;
        } else {
          gx = rad*cos(a);
          gy = rad*sin(a);
          a += adelta;
        }
        vertex(cx+gx,cy+gy);
      }
    
    }
    
    void setup(){
      size(650,650);
      frameRate(24);
    
      smooth();
    
      a = 1;
      b = 2;
      c = 3;
      d = 2;
    
      theta = 0.0;
      dtheta = 8*PI/180;
    
      // declare Input point, Output points
      IP = new float[2];
      OP1 = new float[2];
      OP2 = new float[2];
      // declare ground pivot 1 and 2
      GP1 = new float[2];
      GP2 = new float[2];
    
      if(pdfOutput){
        pdf = (PGraphicsPDF)beginRecord(PDF, "deneme2.pdf");
        beginRecord(pdf);
      }
    
    }
    
    float[] ab(float d1, float d2, float d){
      /*
      finds normalized coordinates (alpha,beta) of point which is at a distance
      d1,d2 from  two points that are located with a distance d apart from each other.
      */
      float[] result = new float[2];
      float delta = (d1-d2)/d;
      float sigma = (d1+d2)/d;
    
      result[0] = (1+sigma*delta)/2.0;
      result[1] = 0.5*sqrt((1-sq(delta))*(sq(sigma)-1));
    
      return result;
    }
    
    void drawConnectingRodx(float r1, float r2, float R1, float R2, float rho1, float rho2, float d, boolean filled){
      /*
      draws connecting rod schematically aligned with x axis.
      (R1,r1) and (R2,r2) are big and small radii at the two ends.
      d is the distance between centers.
      require:
        R1>r1>0,   R2>r2>0,   rho1>0,   rho2>0
      */
      float a = asin((r2-r1)/d);
      float p1x = -r1*sin(a);
      float p1y = +r1*cos(a);
    
      float p2x = d-r2*sin(a);
      float p2y =  +r2*cos(a);
      strokeWeight(0.1);
      point(0.0,0.0);
      point(d,0.0);
    
      float q1x = p1x + sqrt(sq(R1)-sq(r1))*cos(a);
      float q1y = p1y + sqrt(sq(R1)-sq(r1))*sin(a);
    
      float q2x = p2x - sqrt(sq(R2)-sq(r2))*cos(a);
      float q2y = p2y - sqrt(sq(R2)-sq(r2))*sin(a);
    
      float p1s1 = sqrt( sq(R1+rho1)-sq(rho1+r1) );
      float p2s2 = sqrt( sq(R2+rho2)-sq(rho2+r2) );
    
      float b1 = acos((rho1+r1)/(rho1+R1));
      float b2 = acos((rho2+r2)/(rho2+R2));
    
      float w1x = (rho1+R1)*cos(PI/2+a-b1);
      float w1y = (rho1+R1)*sin(PI/2+a-b1);
    
      float w2x = d + (rho2+R2)*cos(PI/2+a+b2);
      float w2y =   + (rho2+R2)*sin(PI/2+a+b2);
    
      float s1x = p1x + p1s1*cos(a);
      float s1y = p1y + p1s1*sin(a);
    
      float s2x = p2x - p2s2*cos(a);
      float s2y = p2y - p2s2*sin(a);
    
      if(filled){
        float dela = 0.03;
    
        fill(100,100,200,150);
        beginShape();
        circularArc(0.,0.,R1,2*PI-(PI/2+a-b1),PI/2+a-b1,dela);
        circularArc(w1x,w1y,rho1,3*PI/2+a-b1,3*PI/2+a,dela);
    
        circularArc(w2x,w2y,rho2,3*PI/2+a,3*PI/2+a+b2,dela);
        circularArc(d,0,R2,PI/2+a+b2,-(PI/2+a+b2),dela);
    
        circularArc(w2x,-w2y,rho2,-(3*PI/2+a+b2),-(3*PI/2+a),dela);
        endShape(CLOSE);
      } else {
        noFill();
        strokeWeight(0.015);
        arc(0,0,        2*R1,2*R1,        PI/2+a-b1,PI); // red-left-top
        arc(w1x,w1y,    2*rho1,2*rho1,    3*PI/2+a-b1,3*PI/2+a); // blue-left-top
    
        line(s1x,s1y,s2x,s2y); // green-top
    
        arc(w2x,w2y,    2*rho2,2*rho2,   3*PI/2+a,3*PI/2+a+b2); //blue-right-top
        arc(d,0,        2*R2,2*R2,       0,PI/2+a+b2); // red-right-top
    
        arc(d,0,        2*R2,2*R2,       -(PI/2+a+b2),0); // red-right-below
        arc(w2x,-w2y,   2*rho2,2*rho2,   -a+PI/2-b2,-a+PI/2);  // blue-right-below
    
        arc(w1x,-w1y,   2*rho1,2*rho1,    PI/2-a,PI/2-a+b1); // blue-left-below
        arc(0,0,        2*R1,2*R1,        PI,3*PI/2+b1-a); // red-left-below
    
        line(s1x,-s1y,s2x,-s2y); // green-below
    
        strokeWeight(0.005);
        line(0.0, 0.0, d, 0.0);
      }
    }
    
    void drawConnectingRod(float r1, float r2, float R1, float R2, float rho1, float rho2, float p1[], float p2[], boolean filled){
      float dx = p2[0]-p1[0];
      float dy = p2[1]-p1[1];
      float angle = atan2(dy,dx);
      float d = sqrt(dx*dx+dy*dy);
    
      pushMatrix();
      translate(p1[0],p1[1]);
      rotate(angle);
      drawConnectingRodx( r1, r2, R1, R2, rho1, rho2, d, filled);
      popMatrix();
    }
    
    void draw(){
    
      translate(width/2, height/2);
      scale(width/2/2, -height/2/2);
    
      background(120,120,200);
      strokeWeight(0.05);
      strokeCap(ROUND);
      strokeJoin(ROUND);
    
      stroke(255,120,127);
    
      // increment current input crank angle
      theta += dtheta;
    
      // set input pivot~point
      IP[0] = a*cos(theta);
      IP[1] = a*sin(theta);
    
      // set output ground pivot~point (constant)
      GP2[0] = d;
      GP2[1] = 0;
    
      // compute the distance between the input pivot and output ground pivot
      float D = dist(IP[0],IP[1],GP2[0],GP2[1]);
    
      // compute local normalized coordinates that describe the point having distances
      // b and c to input pivot and output ground pivot
      float[] alphabeta = ab(b,c,D);
      // compute the two alternates for output pivot
      // alternate 1
      OP1[0] = IP[0] + alphabeta[0]*(GP2[0]-IP[0])+alphabeta[1]*(-(GP2[1]-IP[1]));
      OP1[1] = IP[1] + alphabeta[0]*(GP2[1]-IP[1])+alphabeta[1]*(+(GP2[0]-IP[0]));
      // alternate 2
      OP2[0] = IP[0] + alphabeta[0]*(GP2[0]-IP[0])-alphabeta[1]*(-(GP2[1]-IP[1]));
      OP2[1] = IP[1] + alphabeta[0]*(GP2[1]-IP[1])-alphabeta[1]*(+(GP2[0]-IP[0]));
    
      translate(-0.5,-0.95);
    
      stroke(200,50,50);
      drawConnectingRod(0.1, 0.2, 0.2, 0.3, 0.2, 0.2, GP1, IP,false);
      stroke(50,200,250);
      drawConnectingRod(0.1, 0.2, 0.2, 0.3, 0.2, 0.2, IP , OP1,false);
      stroke(50,50,200);
      drawConnectingRod(0.1, 0.2, 0.2, 0.3, 0.2, 0.2, OP1, GP2,false);
    
      if(pdfOutput){
        pdf.nextPage();
      }
    
    }
    
    void keyPressed(){
      if (key=='q') {
        if(pdfOutput){
          endRecord();
        }
        exit();
      }
    }
    

    One frame from the animation is shown below

    output of "Processing" code for fourbar animation

    References:

    [1] “Einfuhrung in die Getriebelehre” “Hanfried Kerle,Reinhard Pittschellis, Burkhard Corves
    [2] “Theoretical Kinematics” “O. Bottema and B.Roth”

    Some Experiments in Renderman using Aqsis and Pixie

    January 2, 2010

    In this post, simple examples of scenes and shaders in Renderman, that i experimented while trying to learn Renderman, are given. Renderers i experimented are Aqsis and Pixie. First tutorial that i follow is RManNotes (by Steve May) describing how to write shaders.

    RManNotes Tutorial

    I copied two shaders in first section “surf1_1.sl”. For these shaders to work, an include file named “rmannotes.sl” must be also available.

    /* FILE: surf1_1.sl */
    #include "rmannotes.sl"
    
    surface surf1_1()
    {
      color surface_color, layer_color;
      color layer_opac;
    
      /* background layer (layer 0) */
    
      surface_color = color (0.3, 0.3, 0.3); /* grey */
    
      /* layer #1 */
    
      layer_color = color (0,0,1);  /* blue */
      layer_opac = pulse(0.35, 0.65, 0.02, s);
      surface_color = blend(surface_color, layer_color, layer_opac);
    
      /* layer #2 */
    
      layer_color = color (0,1,0);  /* green */
      layer_opac = pulse(0.35, 0.65, 0.0, t);
      layer_opac *= 0.5;
      surface_color = blend(surface_color, layer_color, layer_opac);
    
      /* output */
    
      Ci = surface_color;
    }
    

    To use a shader file, it must be first compiled. In Aqsis and Pixie shader files (*.sl) are compiled by following commands

    rem compiling shader file in Aqsis
    aqsl surf1_1.sl
    rem compiling shader file in Pixie
    sdrc surf1_1.sl
    

    Compilation of shader file (*.sl) by Aqsis and Pixie results in *.slx and *.sdr files, respectively. If shader files are successfully compiled, then they are ready to be used in scene (*.rib) files. To test these shaders, a scene (rib) file containing a sphere and a torus is created.

    #Scene1.rib
    Display "Scene1.png" "file" "rgb"
    Projection "perspective" "fov" 40
    Format 400 400 1
    Translate 0 0 7
    Rotate -120 1 0 0
    #Rotate 60 0 1 0
    WorldBegin
    	LightSource "pointlight" 1 "intensity" 45 "from" [2 -3 4]
    	LightSource "ambientlight" 2 "intensity" 0.5
    	Opacity [0.5 0.5 0.5]
    	Surface "surf1_1"
    	Color 0.5 0.5 1
    	Sphere 0.5 -0.5 0.5 360
    	Torus 1.45 0.35 0 360 360
    WorldEnd
    

    Scene files are rendered by Aqsis and Pixie by following commands

    rem render Scene1.rib by Aqsis
    aqsis Scene1.rib
    rem render Scene1.rib by Pixie
    rndr Scene1.rib
    

    Executing these commands results in png image file “Scene1.png” as indicated by 2nd line of rib file Display "Scene1.png" "file" "rgb".

    The tutorial then gives a variation of above shader in which the stripe is multiplied and diffuse-specular-ambient effects are added

    /* FILE: surf1_2.sl */
    #include "rmannotes.sl"
    
    surface surf1_2(float Ka = 0.7, Kd = 0.7, Ks = 0.3, roughness = 0.1;
    		color specularcolor = 1)
    {
      color surface_color, layer_color;
      color surface_opac, layer_opac;
      float ss, tt;
      point Nf, V;
    
      /* background layer (layer 0) */
    
      surface_color = color (0.3, 0.3, 0.3); /* grey */
      surface_opac = Os;
    
      /* layer #1 */
    
      ss = repeat(s, 4);
      tt = repeat(t, 4);
      layer_color = color (0,0,1);  /* blue */
      layer_opac = pulse(0.35, 0.65, 0.02, ss);
      surface_color = blend(surface_color, layer_color, layer_opac);
      surface_opac = union(surface_opac, layer_opac);
    
      /* layer #2 */
    
      ss = repeat(s, 4);
      tt = repeat(t, 4);
      layer_color = color (0,1,0);  /* green */
      layer_opac = pulse(0.35, 0.65, 0.02, tt);
      layer_opac *= 0.5;
      surface_color = blend(surface_color, layer_color, layer_opac);
      surface_opac = union(surface_opac, layer_opac);
    
      /* illum */
    
      Nf = faceforward(normalize(N), I);
      V = -normalize(I);
      surface_color = surface_color * (Ka * ambient() + Kd * diffuse(Nf))
        + specularcolor * Ks * specular(Nf, V, roughness);
    
      /* output */
    
      Oi = surface_opac;
      Ci = surface_color * surface_opac;
    }
    

    and rib file is modified as follows

    #Scene2.rib
    Display "Scene2.png" "file" "rgb"
    Projection "perspective" "fov" 40
    ...
    ...
    	Surface "surf1_2"
    

    new render with new shader “surf1_2.sl” is obtained via

    aqsl surf1_2.sl
    aqsis Scene2.rib
    

    resulting in following image

    As a digression from the tutorial, at this point, a silhouette shader is attempted. Condition for a surface point to be a silhouette point is that normal of surface is perpendicular to incident ray, both of which available in renderman shading language as global variables. Silhouette condition is expressed as

    \displaystyle \mathbf{V\cdot N} = 0

    Since surface color must be at peak this is reversed by

    \displaystyle 1-\mathbf{V\cdot N}

    and lower values are rejected by utility macro “pulse” from

    /* macro pulse from "rmannotes.sl" */
    #define pulse(a,b,fuzz,x) (smoothstep((a)-(fuzz),(a),(x)) - \
    			   smoothstep((b)-(fuzz),(b),(x)))
    

    in this example rejection range is selected as 0-0.7

    complete shader code for silhouette is

    #include "rmannotes.sl"
    
    surface surf_silhouette()
    {
      point Nf, V;
      color edgecolor = color "rgb" (0.0, 0.0, 0.0);
    
      Nf = faceforward(normalize(N), I);
      V = -normalize(I);
    
      Oi = Os;
      Ci = blend( edgecolor, Cs, pulse(0.7, 1.0, 0.01, 1-abs(V.Nf) ));
    }
    

    rib file is modified as

    #SceneSilhouette.RIB
    Display "SceneSilhouette.png" "file" "rgb"
    ...
    ...
    	Surface "surf_silhouette"
    

    The resulting image is given below

    Another attempt is to use Screen shader by Larry Gritz

    If Aqsis is used on scene (rib) file given above and screen shader, an error message saying ” token density is not defined” appears. This problem is mentioned in companion site for book “Rendering for Beginners” by Saty Raghavachary, saying that all argument names in a call to shader must be preceded by argument types (for many renderers including Aqsis). With this correction new rib file is as follows

    #SceneScreen.RIB
    Display "SceneScreen.png" "file" "rgb"
    Projection "perspective" "fov" 40
    Format 640 640 1
    Translate 0 0 7
    Rotate -120 1 0 0
    WorldBegin
    	LightSource "pointlight" 1 "intensity" 45 "from" [2 -3 4]
    	LightSource "ambientlight" 2 "intensity" 0.75
    	Opacity [0.5 0.5 0.5]
    	Surface "screen_aa" "float Ka" 1 "float Kd" 0.75 "float Ks" 0.4 "float roughness" 0.1
    	 "color specularcolor" [1 1 1] "float density" 0.25 "float frequency" 20
    	Color 0.5 0.5 1 
    	Sphere 0.5 -0.5 0.5 360
    	Color 0.9 0.5 1
    	Torus 1.45 0.35 0 360 360
    WorldEnd
    

    Screen shader code is also given below for completeness (copied from RenderGibbon by Mark Schmelzenbach

    /*
     * screen_aa.sl -- RenderMan compatible shader for a metalic screen.
     *
     * DESCRIPTION:
     * Makes a surface that looks like a metal screen. Strips of metal run
     * parallel to lines of s and t.  You can adjust the Ka, Kd, Ks, etc.
     * to change the material appearance.  This texture antialiases pretty
     * well, even with only one sample per pixel.
     *
     * PARAMETERS:
     * Ka, Kd, Ks, roughness, specularcolor - work like the plastic shader
     *   frequency - how many cycles of screen in st space
     *   density - how much of each cycle is opaque?
     *
     * AUTHOR: written by Larry Gritz
     */ 
    
    #define boxstep(a,b,x) (clamp(((x)-(a))/((b)-(a)),0,1))
    #define MINFILTERWIDTH 1.0e-7
    
    surface screen_aa (float Ka = 1, Kd = 0.75, Ks = 0.4, roughness = 0.1;
       color specularcolor = 1;float density = 0.25, frequency = 20;)
    {
      point Nf;     /* Forward facing Normal vector */
      point IN;     /* normalized incident vector */
      float d;      /* Density at the sample point */
      float ss, tt; /* s,t, parameters in phase */
      float swidth, twidth, GWF, w, h;
    
      /* Compute a forward facing normal */
      IN = normalize (I);
      Nf = faceforward (normalize(N), I);
    
      /* Determine how wide in s-t space one pixel projects to */
      swidth = max (abs(Du(s)*du) + abs(Dv(s)*dv), MINFILTERWIDTH) *
    frequency;
      twidth = max (abs(Du(t)*du) + abs(Dv(t)*dv), MINFILTERWIDTH) *
    frequency;
    
      /* Figure out where in the pattern we are */
      ss = mod (frequency * s, 1);
      tt = mod (frequency * t, 1);
    
      /* Figure out where the strips are. Do some simple antialiasing. */
      GWF = density*0.5;
      if (swidth >= 1)
          w = 1 - 2*GWF;
      else w = clamp (boxstep(GWF-swidth,GWF,ss), max(1-GWF/swidth,0), 1)
     - clamp (boxstep(1-GWF-swidth,1-GWF,ss), 0, 2*GWF/swidth);
      if (twidth >= 1)
          h = 1 - 2*GWF;
      else h = clamp (boxstep(GWF-twidth,GWF,tt), max(1-GWF/twidth,0),1)
     - clamp (boxstep(1-GWF-twidth,1-GWF,tt), 0, 2*GWF/twidth);
      /* This would be the non-antialiased version:
       *    w = step (GWF,ss) - step(1-GWF,ss);
       *    h = step (GWF,tt) - step(1-GWF,tt);
       */
      d = 1 - w*h;
    
      Oi = d;
      if (d > 0) {
          Ci = Oi * ( Cs * (Ka*ambient() + Kd*diffuse(Nf)) +
     specularcolor * Ks*specular(Nf,-IN,roughness));
        }
      else
          Ci = 0;
    }
    

    Result of this shader on Scene1.rib is given below

    Using TAPENADE

    November 27, 2009

    TAPENADE is a Automatic Differentiation (AD) tool for FORTRAN. It generates differentiation code by source code transformation. This post gives examples of TAPENADE usage.

    TAPENADE Example

    For illustration of TAPENADE usage, following test expression is selected

    \displaystyle z\left(x,y\right)=\sin \left( x+2y\right) \cos \left( x-\sin \left( y\right) \right)

    This expression/function is coded in FORTRAN as a subroutine which takes 3 arguments whose in/out behaviour is stated by intent directives

  • x,y : input arguments
  • z : output arguments
  • ! FILE : test1.f90
    subroutine test1(x,y,z)
    	real, intent(in) :: x,y
    	real, intent(out) :: z
    
    	z = sin(x+2*y)*cos(x-sin(y));
    end subroutine test1
    

    TAPENADE executed to differentiate test1 subroutine (in tangent mode) generates following differentiation code

    !        Generated by TAPENADE     (INRIA, Tropics team)
    !  Tapenade 3.3 (r3163) - 09/25/2009 09:03
    !
    !  Differentiation of test1 in forward (tangent) mode:
    !   variations  of output variables: z
    !   with respect to input variables: x y
    SUBROUTINE TEST1_D(x, xd, y, yd, z, zd)
      IMPLICIT NONE
      REAL, INTENT(IN) :: x, y
      REAL, INTENT(IN) :: xd, yd
      REAL, INTENT(OUT) :: z
      REAL, INTENT(OUT) :: zd
      REAL :: arg1
      REAL :: arg1d
      INTRINSIC COS
      INTRINSIC SIN
      arg1d = xd - yd*COS(y)
      arg1 = x - SIN(y)
      zd = (xd+2*yd)*COS(x+2*y)*COS(arg1) - SIN(x+2*y)*arg1d*SIN(arg1)
      z = SIN(x+2*y)*COS(arg1)
    END SUBROUTINE TEST1_D
    

    Note that each input and output variable that user specified caused insertion of “differentials” of those variables, that is, the input code

    subroutine test1(x,y,z)

    results in

    SUBROUTINE TEST1_D(x, xd, y, yd, z, zd)

    In this example,

  • differentials dx,dy for input arguments x,y
  • differential dz for output argument z
  • are inserted just after the original argument.

    Actual differential of test1 function is (total differential)

    \displaystyle dz=\frac{\partial z}{\partial x}dx+\frac{\partial z}{\partial y}dy

    with

    \frac{\partial z}{\partial x}=\cos \left( x+2y\right) \cos \left( x-\sin \left( y\right) \right) -\sin \left( x+2y\right) \sin \left( x-\sin \left( y\right) \right)
    \frac{\partial z}{\partial y}=2\cos \left( x+2y\right) \cos \left( x-\sin \left( y\right) \right) +\sin \left( x+2y\right) \sin \left( x-\sin \left( y\right) \right) \cos \left( y\right)

    Note that in order to compute partial derivatives \frac{\partial z}{\partial x} and \frac{\partial z}{\partial y}, the generated derivative routine test1_d must be called twice

    ! derivative computation point
    x = ...
    y = ...
    xd = 1.0
    yd = 0.0
    CALL TEST1_D(x, xd, y, yd, z, zd) ! z and zd ~ dz/dx is obtained
    xd = 0.0
    yd = 1.0
    CALL TEST1_D(x, xd, y, yd, z, zd) ! z and zd ~ dz/dy is obtained
    

    In order to avoid such multiple calls, TAPENADE has an alternate differentiation mode called “Tangent Vectorial Mode”

    If TAPENADE code is generated with this option then following differentiation code is obtained

    !        Generated by TAPENADE     (INRIA, Tropics team)
    !  Tapenade 3.3 (r3163) - 09/25/2009 09:03
    !
    !  Differentiation of test1 in forward (tangent) mode: (multi-directional mode)
    !   variations  of output variables: z
    !   with respect to input variables: x y
    SUBROUTINE TEST1_DV(x, xd, y, yd, z, zd, nbdirs)
      USE DIFFSIZES
    !  Hint: nbdirsmax should be the maximum number of differentiation directions
      IMPLICIT NONE
      REAL, INTENT(IN) :: x, y
      REAL, INTENT(IN) :: xd(nbdirsmax), yd(nbdirsmax)
      REAL, INTENT(OUT) :: z
      REAL, INTENT(OUT) :: zd(nbdirsmax)
      REAL :: arg1
      REAL :: arg1d(nbdirsmax)
      REAL :: arg2
      REAL :: arg2d(nbdirsmax)
      INTEGER :: nd
      INTEGER :: nbdirs
      INTRINSIC COS
      INTRINSIC SIN
      arg1 = x + 2*y
      arg2 = x - SIN(y)
      DO nd=1,nbdirs
        arg1d(nd) = xd(nd) + 2*yd(nd)
        arg2d(nd) = xd(nd) - yd(nd)*COS(y)
        zd(nd) = arg1d(nd)*COS(arg1)*COS(arg2) - SIN(arg1)*arg2d(nd)*SIN(&
    &      arg2)
      END DO
      z = SIN(arg1)*COS(arg2)
    END SUBROUTINE TEST1_DV
    

    Note that generated code executes differentiation code for all given differentiation directions inside a loop. This, in effect, corresponds to calling “single-direction” mode successively, but with the direction loop inside the differentiation code, subroutine arguments are not passed to calls again and again, increasing the efficiency of code. Another benefit of this mode is hiding details of successive calls.

    Following (incomplete) program illustrates use of “multi-directional” mode.

    integer, parameter :: nbdirs = 8
    real*8 xd(nbdirs),yd(nbdirs)
    xd = (/1.0, 0.0/)
    yd = (/0.0, 1.0/)
    call TEST1_DV(x, xd, y, yd, z, zd, nbdirs) ! zd contains all partial derivatives
    

    TAPENADE GUI Usage

    First click add button

    Select the source fi le(s) test.f90

    In listbox component all source code added is shown

    Click Parse Button. On the left drop-down list box, subroutines/functions available in parsed source codes will appear

    In input variables, click Add and select x,y

    Pressing OK will add selected variable(s) to input variables listbox component

    In output variable, click Add and select z

    Pressing OK will add selected variable(s) to output variables listbox component

    Finally select differentiation mode and click differentiate. This will generate code and automatically open html report page (if browser is properly defi ned in environment)

    Reverse Mode

    Reverse mode first illustrated by a sample routine from INRIA report RR-5237 . This routine implements function

    \displaystyle o_1 = z\left(  x^{2}+1\right)  y^{2} if \displaystyle \left(  x^{2}+1\right)  y^{2} < 10
    \displaystyle o_1 = 10z otherwise

    Exact gradient of this expression is

    \displaystyle \frac{\partial o_{1}}{\partial\left( x,y,z\right) } = \left[ 2zxy^{2}, z\left(  x^{2}+1\right), \left(  x^{2}+1\right)  y^{2}\right]

    if \displaystyle \left(  x^{2}+1\right)  y^{2} < 10

    subroutine sub1(x,y,z,o1)
    	real, intent(in) :: x,y,z
    	real, intent(out) :: o1
    
    	x = y*x
    	r = x*x+y*y   ! x**2*y**2+y**2
    	if (r>10.0) then
    		r = 10.0
    	end if
    
    	o1 = r*z      ! z*(x**2*y**2+y**2)
    end subroutine sub1
    

    differentiating this routine in reverse mode using TAPENADE yields

    !        Generated by TAPENADE     (INRIA, Tropics team)
    !  Tapenade 3.3 (r3163) - 09/25/2009 09:03
    !
    !  Differentiation of sub1 in reverse (adjoint) mode:
    !   gradient, with respect to input variables: x y z o1
    !   of linear combination of output variables: o1
    SUBROUTINE SUB1_B(x, xb, y, yb, z, zb, o1, o1b)
      IMPLICIT NONE
      REAL :: x, y, z
      REAL :: xb, yb, zb
      REAL :: o1
      REAL :: o1b
      REAL :: r
      REAL :: rb
      INTEGER :: branch
      CALL PUSHREAL4(x)
      x = y*x
      r = x*x + y*y
      IF (r .GT. 10.0) THEN
        r = 10.0
        CALL PUSHINTEGER4(1)
      ELSE
        CALL PUSHINTEGER4(0)
      END IF
      rb = z*o1b
      zb = r*o1b
      CALL POPINTEGER4(branch)
      IF (.NOT.branch .LT. 1) rb = 0.0
      xb = 2*x*rb
      CALL POPREAL4(x)
      yb = x*xb + 2*y*rb
      xb = y*xb
      o1b = 0.0
    END SUBROUTINE SUB1_B
    

    Note that TAPENADE does not include evaluation of original function in reverse-mode. Some literature states that this TAPENADE specific.

    Driver program for reverse differentiation code is given below. Program computes gradient at x=1,y=1.5,z=2 which bypass the if-statement in the code. The gradient of linear combination selected output variables o1 is returned in differentials of input variables xb,yb,zb. Weightings are specified in differentials of output variables o1b. So since in this example there is one output variable, o1b=1 yields gradient of o1 in xb,yb,zb.

    program usesub1b
    	implicit none
    	real x,xb,y,yb,z,zb,o1,o1b
    	o1b = 1.0
    	x = 1.0
    	y = 1.5
    	z = 2.0
    	call SUB1_B(x, xb, y, yb, z, zb, o1, o1b)
    
    	print *, "TAPENADE computed gradient:"
    	write(6,'(3(F12.6))') xb,yb,zb
    
    	print *, "exacy gradient (without if):"
    	write(6,'(3(F12.6))') z*2*x*y**2, z*(x**2+1)*2*y, (x**2+1)*y**2
    end program usesub1b
    

    Note that reverse mode differentiation code contains PUSH/POP (store value-retrieve value) operations. To my knowledge, this is the case when some branching like if-statement, subroutine-call occurs in “source” code. Again to my knowledge, number of occurrences of PUSH/POP operations is dependent on the AD tool reverse mode differentiation strategies “store-all”, “re-compute-all”,… . These PUSH/POP operations are defined on separate source files (ADFirstAidKit.tar.gz) available for download at TAPENADE site. Compilation is described in TAPENADE tutorial. For this example, only adBuffer.f,adStack.c is needed. Following the instructions in tutorial code is compiled by following commands

    rem from http://www-sop.inria.fr/tropics/tapenade/tutorial.html
    gfortran -c usesub1b.f90
    gfortran -c sub1_b.f90
    gfortran -c adBuffer.f
    gcc -c adStack.c
    
    gfortran sub1_b.o adBuffer.o adStack.o usesub1b.o -o usesub1b
    

    output of executable is as follows

     TAPENADE computed gradient:
        9.000000   12.000000    4.500000
     exacy gradient (without if):
        9.000000   12.000000    4.500000
    

    Second example is given to illustrate concept of gradient of linear combination of output variables. In this example two input variables map to \left(x,y\right) \rightarrow \left(\begin{array}[c]{c}o_{1}\\o_{2}\end{array}\right)  =\left( \begin{array} [c]{c} \frac{x^{2}-y^{2}}{2}\\ \frac{x^{2}+y^{2}}{2} \end{array}\right) . This function code is given below

    subroutine samplefun(x,y,z1,z2)
    	real, intent(in) :: x,y
    	real, intent(out) :: z1,z2
    
    	z1 = 0.5*(x**2-y**2)
    	z2 = 0.5*(x**2+y**2)
    
    end subroutine samplefun
    

    and its reverse-mode differentiation is obtained as

    !        Generated by TAPENADE     (INRIA, Tropics team)
    !  Tapenade 3.3 (r3163) - 09/25/2009 09:03
    !
    !  Differentiation of samplefun in reverse (adjoint) mode:
    !   gradient, with respect to input variables: x y z1 z2
    !   of linear combination of output variables: z1 z2
    SUBROUTINE SAMPLEFUN_B(x, xb, y, yb, z1, z1b, z2, z2b)
      IMPLICIT NONE
      REAL, INTENT(IN) :: x, y
      REAL :: xb, yb
      REAL :: z1, z2
      REAL :: z1b, z2b
      xb = 0.5*2*x*z1b + 0.5*2*x*z2b
      yb = 0.5*2*y*z2b - 0.5*2*y*z1b
      z1b = 0.0
      z2b = 0.0
    END SUBROUTINE SAMPLEFUN_B
    

    Define linear combination of two outputs o_1,o_2 as J=w_{1}o_{1}+w_{2}o_{2}=w_{1}\frac{x^{2}-y^{2}}{2}+w_{2}\frac{x^{2}+y^{2}}{2} then gradient of this new quantity is

    \nabla J=\left[ \begin{array} [c]{c} \left(  w_{1}+w_{2}\right)  x\\ \left(  w_{2}-w_{1}\right)  y \end{array} \right]

    To check this in differentiation code, manually call SAMPLEFUN_B with z1b=w1, z2b=w2 to get

    xb = 0.5*2*x*z1b + 0.5*2*x*z2b => 0.5*2*x*w1 + 0.5*2*x*w2 \rightarrow x w_1+x w_2
    yb = 0.5*2*y*z2b - 0.5*2*y*z1b => 0.5*2*y*z2b - 0.5*2*y*z1b \rightarrow y w_2-y w_1

    which is recognized as actual gradient as required.

    Missing functions

    TAPENADE handles all common mathematical functions sin, cos, exp, sqrt in the differentiation process, but if some function is not known to TAPENADE and/or its implementation is missing, instead of giving an error and quitting, TAPENADE continues with generation of differentiation code using the apparent signature of the function. Then user is responsible for supplying differentiation code of the function with missing implementation. To exemplify situation, the function

    \displaystyle s\left(\boldsymbol{r}\right)=\frac{\sin\left(\sqrt{\boldsymbol{r}\cdot\boldsymbol{r}}\right)}{\sqrt{\boldsymbol{r}\cdot\boldsymbol{r}}}

    that computes the sinc-of-radial distance in nD space (using the DOT_PRODUCT) is differentiated with TAPENADE. Code to be differentiated is given below

    ! FILE:radsinc.f90
    subroutine radsinc(r,n,rs)
    	integer, intent(in) :: n
    	real, intent(in) :: r(n)
    	real, intent(out) :: rs
    
    	real rad
    
    	rad = sqrt(dot_product(r,r))
    	if (rad==0) then
    		rs = 1.0
    	else
    		rs = sin( rad ) / rad
    	end if
    end subroutine radsinc
    

    and differentiation code in tangent mode is obtained as

    ! FILE: radsinc_d.f90
    !        Generated by TAPENADE     (INRIA, Tropics team)
    !  Tapenade 3.3 (r3163) - 09/25/2009 09:03
    !
    !  Differentiation of radsinc in forward (tangent) mode:
    !   variations  of output variables: r rs
    !   with respect to input variables: r
    SUBROUTINE RADSINC_D(r, rd, n, rs, rsd)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: n
      REAL, INTENT(IN) :: r(n)
      REAL, INTENT(IN) :: rd(n)
      REAL, INTENT(OUT) :: rs
      REAL, INTENT(OUT) :: rsd
      REAL :: rad
      REAL :: radd
      EXTERNAL DOT_PRODUCT
      EXTERNAL DOT_PRODUCT_D
      REAL :: DOT_PRODUCT
      REAL :: DOT_PRODUCT_D
      REAL :: result1
      REAL :: result1d
      INTRINSIC SIN
      INTRINSIC SQRT
      result1d = DOT_PRODUCT_D(r, rd, r, rd, result1)
      IF (result1 .EQ. 0.0) THEN
        radd = 0.0
      ELSE
        radd = result1d/(2.0*SQRT(result1))
      END IF
      rad = SQRT(result1)
      IF (rad .EQ. 0) THEN
        rs = 1.0
        rsd = 0.0
      ELSE
        rsd = (radd*COS(rad)*rad-SIN(rad)*radd)/rad**2
        rs = SIN(rad)/rad
      END IF
    END SUBROUTINE RADSINC_D
    

    Actual differential of the function is

    \displaystyle ds =\frac{\boldsymbol{r}\cdot\boldsymbol{dr}}{\left(\boldsymbol{r}\cdot\boldsymbol{r}\right)^{\frac{3}{2}}}\left(\cos\left(\sqrt{\boldsymbol{r}\cdot\boldsymbol{r}}\right)\sqrt{\boldsymbol{r}\cdot\boldsymbol{r}}-\sin\left(\sqrt{\boldsymbol{r}\cdot\boldsymbol{r}}\right)\right)

    Since the function DOT_PRODUCT is not in TAPENADE standard function list and no implementation is available, TAPENADE cannot generate a differentiation code for function DOT_PRODUCT. So it just emits the differential code result1d = DOT_PRODUCT_D(r, rd, r, rd, result1) (line 25) and expects the user to supply this function at compile time. For this reason, it informs the compiler that actual object code is defined elsewhere by inserting EXTERNAL declaration(s) in the source code (line 18). Derivative of dot product can be obtained from

    \displaystyle \frac{d}{dt}\left(\boldsymbol{a\cdot b}\right)=\frac{d\boldsymbol{a}}{dt}\cdot\boldsymbol{b}+\boldsymbol{a}\cdot\frac{d\boldsymbol{b}}{dt}

    which is translated to code as follows

    ! FILE: DOT_PRODUCT_D.f90
    real function dot_product_d(a, ad, b, bd)
    	implicit none
    	real, intent(in), dimension(:) :: a, b, ad, bd
    	dot_product_d = dot_product(ad,b) + dot_product(a,bd)
    end function dot_product_d
    

    From the generated code, it is observed that TAPENADE generated code computes function and its derivative in the same routine, so if conformance to originally generated code is sought then above definition is updated as

    ! FILE: DOT_PRODUCT_D.f90
    real function dot_product_d(a, ad, b, bd,dotprod)
    	implicit none
    	real, intent(in), dimension(:) :: a, b, ad, bd
    	real, intent(out) :: dotprod
    	dotprod = dot_product(a,b)
    	dot_product_d = dot_product(ad,b) + dot_product(a,bd)
    end function dot_product_d
    

    To test the TAPENADE generated routine following test program is used

    ! FILE: main.f90
    program main
    	implicit none
    
    	interface
    		real function dot_product_d(a, ad, b, bd,dotprod)
    			implicit none
    			real, intent(in), dimension(:) :: a, b, ad, bd
    			real, intent(out) :: dotprod
    		end function dot_product_d
    	end interface
    
    	real, dimension(4) :: r, rd
    	real :: rs, rsd, rstrue, rsdtrue, rdotr, rdotrd, rdotr2
    	real ::  dotprod,dotprod_d
    
    	r = (/ 1., 2., 3., 4. /)		! vector in 4D
    	rd = (/ 40, 30, 20, 10 /)	! rate(/differential) of vector
    
    	! check dot product
    	rdotr = dot_product(r,r)
    	rdotrd = dot_product(r,rd)
    	print *, "r-dot-dr:", rdotrd
    	print *, "r-dot-r :", rdotr
    
    	! check derivative of dot product
    	dotprod_d = dot_product_d(r, rd, r, rd, rdotr2)
    	print *, "r-dot-r   :", rdotr2
    	print *, "d(r-dot-r):", dotprod_d
    
    	! call TAPENADE generated routine
    	! SUBROUTINE RADSINC_D(r, rd, n, rs, rsd)
    	call RADSINC_D(r, rd, 4, rs, rsd)
    
    	! compute analytical result
    	rsdtrue = (rdotrd/rdotr**1.5)*(cos(sqrt(rdotr))*sqrt(rdotr)-sin(sqrt(rdotr)))
    	rstrue = sin(sqrt(rdotr))/sqrt(rdotr)
    
    	print *, "analytical results:"
    	write (*,'(A8,F16.8,X,A8,F16.8)')  "sinc = ",rstrue, "dsinc = ",rsdtrue
    
    	print *, "TAPENADE results:"
    	write (*,'(A8,F16.8,X,A8,F16.8)')  "sinc = ",rs, "dsinc = ",rsd
    end program main
    

    Compilation command is

    gfortran dotproduct_d.f90 radsinc_d.f90 main.f90 -o main.exe
    

    Program tests the routine for the inputs

    r=\left[1,2,3,4\right]\\dr=\left[40,30,20,10\right]

    output of program is given below

     r-dot-dr:   200.00000
     r-dot-r :   30.000000
     r-dot-r   :   30.000000
     d(r-dot-r):   400.00000
     analytical results:
     sinc =      -0.13172643 dsinc =       5.49430418
     TAPENADE results:
     sinc =      -0.13172643 dsinc =       5.49430418
    

    showing that the results agree

    Note that main program contains an interface declaration for dot_product_d. This is required for routines with assumed shaped array arguments as explained in Fortran Bugs/Surprises. Since TAPENADE generated routine also uses the routine dot_product_d same interface is added to SUBROUTINE RADSINC_D .This eliminates some of the previous declarations added by TAPENADE. Modified version of TAPENADE code is given below

    ! FILE: RAD_SINC.f90
    !        Generated by TAPENADE     (INRIA, Tropics team)
    !  Tapenade 3.3 (r3163) - 09/25/2009 09:03
    !
    !  Differentiation of radsinc in forward (tangent) mode:
    !   variations  of output variables: r rs
    !   with respect to input variables: r
    SUBROUTINE RADSINC_D(r, rd, n, rs, rsd)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: n
      REAL, INTENT(IN) :: r(n)
      REAL, INTENT(IN) :: rd(n)
      REAL, INTENT(OUT) :: rs
      REAL, INTENT(OUT) :: rsd
      REAL :: rad
      REAL :: radd
    
      interface
      	real function dot_product_d(a, ad, b, bd,dotprod)
      		implicit none
      		real, intent(in), dimension(:) :: a, b, ad, bd
      		real, intent(out) :: dotprod
    		end function dot_product_d
    	end interface
    
    ! Following declaration no more needed due to above interface definition
    !EXTERNAL DOT_PRODUCT
    !EXTERNAL DOT_PRODUCT_D
    !REAL :: DOT_PRODUCT
    !REAL :: DOT_PRODUCT_D
    
      REAL :: result1
      REAL :: result1d
      INTRINSIC SIN
      INTRINSIC SQRT
      result1d = DOT_PRODUCT_D(r, rd, r, rd, result1)
      IF (result1 .EQ. 0.0) THEN
        radd = 0.0
      ELSE
        radd = result1d/(2.0*SQRT(result1))
      END IF
      rad = SQRT(result1)
      IF (rad .EQ. 0) THEN
        rs = 1.0
        rsd = 0.0
      ELSE
        rsd = (radd*COS(rad)*rad-SIN(rad)*radd)/rad**2
        rs = SIN(rad)/rad
      END IF
    END SUBROUTINE RADSINC_D
    

    Lines 18-24 are the modifications to original code and 27-30 are the commented lines.

    Similar Technique

    For some cases, need for derivative of missing function may be avoided by some minor modifications to source code. For example suppose that there is a routine that assembles a linear system (a matrix and a RHS vector) and solves linear system using a subroutine from compiled library whose source code is not available

    subroutine SOLVE(input1, input2, x, n)
    	integer, intent(in) :: n
    	real, intent(in) :: input1, input2
    	real, intent(out) :: x(n)
    
    	real :: A(n,n), b(n)
    
    	call assemble(input1, input2, A, b, n)
    	call solveLinearSystem(A,b,x)
    end subroutine SOLVE
    

    when differentiated in tangent mode, following code is obtained from TAPENADE

    !        Generated by TAPENADE     (INRIA, Tropics team)
    !  Tapenade 3.3 (r3163) - 09/25/2009 09:03
    !
    !  Differentiation of solve in forward (tangent) mode:
    !   variations  of output variables: x input1 input2
    !   with respect to input variables: input1 input2
    SUBROUTINE SOLVE_D(input1, input1d, input2, input2d, x, xd, n)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: n
      REAL, INTENT(IN) :: input1, input2
      REAL, INTENT(IN) :: input1d, input2d
      REAL, INTENT(OUT) :: x(n)
      REAL, INTENT(OUT) :: xd(n)
      REAL :: a(n, n), b(n)
      REAL :: ad(n, n), bd(n)
      EXTERNAL ASSEMBLE
      EXTERNAL ASSEMBLE_D
      EXTERNAL SOLVELINEARSYSTEM
      EXTERNAL SOLVELINEARSYSTEM_D
      ad = 0.0
      bd = 0.0
      CALL ASSEMBLE_D(input1, input1d, input2, input2d, a, ad, b, bd, n)
      xd = 0.0
      CALL SOLVELINEARSYSTEM_D(a, ad, b, bd, x, xd)
    END SUBROUTINE SOLVE_D
    

    but as stated in beginning, there is no way to obtain or too difficult/impractical to obtain differentiation code SOLVELINEARSYSTEM_D. Trick is to obtain derivatives by some manipulation on the total differential of original expression \mathbf{Ar}=\mathbf{b} \Rightarrow d\mathbf{A}\cdot\mathbf{r}+\mathbf{A}\cdot d\mathbf{r}=d\mathbf{b} which yields

    \displaystyle d\mathbf{r}=\mathbf{A}^{-1}\cdot\left(d\mathbf{b}-d\mathbf{A}\cdot\mathbf{r}\right)

    To utilize this formula, TAPENADE code is modified as follows

    !        Generated by TAPENADE     (INRIA, Tropics team)
    !  Tapenade 3.3 (r3163) - 09/25/2009 09:03
    !
    !  Differentiation of solve in forward (tangent) mode:
    !   variations  of output variables: x input1 input2
    !   with respect to input variables: input1 input2
    SUBROUTINE SOLVE_D(input1, input1d, input2, input2d, x, xd, n)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: n
      REAL, INTENT(IN) :: input1, input2
      REAL, INTENT(IN) :: input1d, input2d
      REAL, INTENT(OUT) :: x(n)
      REAL, INTENT(OUT) :: xd(n)
      REAL :: a(n, n), b(n)
      REAL :: ad(n, n), bd(n)
      EXTERNAL ASSEMBLE
      EXTERNAL ASSEMBLE_D
      EXTERNAL SOLVELINEARSYSTEM
      EXTERNAL SOLVELINEARSYSTEM_D
      ad = 0.0
      bd = 0.0
      CALL ASSEMBLE_D(input1, input1d, input2, input2d, a, ad, b, bd, n)
      xd = 0.0
      !CALL SOLVELINEARSYSTEM_D(a, ad, b, bd, x, xd)
      CALL SOLVELINEARSYSTEM(a, b, x)
      CALL SOLVELINEARSYSTEM(a, bd-matmul(ad,x), xd)
    
    END SUBROUTINE SOLVE_D
    

    Case of MATMUL function

    MATMUL function is similar to DOT_PRODUCT in that both are inner-products, essentially. Former one is acting on second and first order tensors and latter one acting on first order tensors (vectors). So similar to DOT_PRODUCT example above, manual declaration for differentiation of MATMUL would be,

    \displaystyle d\left(\mathbf{A}\cdot\mathbf{b}\right)=d\mathbf{A}\cdot\mathbf{b}+\mathbf{A}\cdot d\mathbf{b}

    which is coded as

    function matmul_d(a,ad,b,bd,ab)
    	implicit none
    	real, intent(in), dimension(:,:) :: a,b,ad,bd
    	real, intent(out), dimension( size(a,1), size(b,2) ) :: ab
    
    	real matmul_d( size(a,1), size(b,2) )
    
    	ab = matmul(a,b)
    	matmul_d = matmul(ad,b) + matmul(a,bd)
    end function matmul_d
    

    Such declaration poses some Fortran 90 problems. In fortran(90) MATMUL can be used to multiply vector-matrix, matrix-matrix and matrix-vector, but above declaration gives signature for matrix-matrix version. One possible solution is to use generic-interfaces but i did not try this yet; instead i devised example consisting of only matrix-matrix multiplications. Example i choose is quadratic form given as

    \displaystyle Q = \mathbf{x}^{\top}\cdot\mathbf{A}\cdot\mathbf{x}

    whose (exact) differential is given by

    \displaystyle dQ = d\mathbf{x}^{\top}\cdot\mathbf{A}\cdot\mathbf{x}+\mathbf{x}^{\top}\cdot d\mathbf{A}\cdot\mathbf{x}+\mathbf{x}^{\top}\cdot\mathbf{A}\cdot d\mathbf{x}

    quadratic function is coded in FORTRAN as follows

    subroutine quadraticForm(A,x,QF)
    	implicit none
    	real, intent(in), dimension(2,2) :: A
    	real, intent(in), dimension(2,1) :: x
    	real, intent(out), dimension(1,1) :: QF
    
    	real :: xp(2,1), xpt(1,2)
    
    	xp = matmul(A,x)
    	xpt = transpose(xp)
    	QF = matmul(xpt,x)
    end subroutine quadraticForm
    

    Note that vectors are replaced with column or row vector versions to conform to MATMUL_D declaration. Code obtained from tangent mode differentiation of TAPENADE is modified to yield following

    !        Generated by TAPENADE     (INRIA, Tropics team)
    !  Tapenade 3.3 (r3163) - 09/25/2009 09:03
    !
    !  Differentiation of quadraticform in forward (tangent) mode:
    !   variations  of output variables: qf x a
    !   with respect to input variables: x a
    SUBROUTINE QUADRATICFORM_D(a, ad, x, xd, qf, qfd)
      IMPLICIT NONE
    
    	INTERFACE
    	 	function matmul_d(a,ad,b,bd,ab)
    			implicit none
    			real, intent(in), dimension(:,:) :: a,b,ad,bd
    			real, intent(out), dimension( size(a,1), size(b,2) ) :: ab
    
    			real matmul_d( size(a,1), size(b,2) )
    		end function matmul_d
    
    		function dot_product_d(a,ad,b,bd,ab)
    			implicit none
    			real dot_product_d
    			real, intent(in), dimension(:) ::  a,ad,b,bd
    			real, intent(out) :: ab
    		end function dot_product_d
    
    		function transpose_d(v, vd, vt)
    			implicit none
    			real, dimension(:,:), intent(in) :: v, vd
    			real, dimension(size(v,2), size(v,1) ), intent(out) :: vt
    			real, dimension( size(v,2), size(v,1) ) :: transpose_d
    
    		end function transpose_d
    
    	END INTERFACE
    
      REAL, DIMENSION(2, 2), INTENT(IN) :: a
      REAL, DIMENSION(2, 2), INTENT(IN) :: ad
      REAL, DIMENSION(2, 1), INTENT(IN) :: x
      REAL, DIMENSION(2, 1), INTENT(IN) :: xd
      REAL, DIMENSION(1, 1), INTENT(OUT) :: qf
      REAL, DIMENSION(1, 1), INTENT(OUT) :: qfd
      REAL :: xp(2, 1), xpt(1, 2)
      REAL :: xpd(2, 1), xptd(1, 2)
      ! EXTERNAL MATMUL
      ! EXTERNAL MATMUL_D
      ! EXTERNAL TRANSPOSE
      ! EXTERNAL TRANSPOSE_D
      ! REAL :: MATMUL
      ! REAL :: MATMUL_D
      ! REAL :: TRANSPOSE
      ! REAL :: TRANSPOSE_D
      xpd = MATMUL_D(a, ad, x, xd, xp)
      xptd = TRANSPOSE_D(xp, xpd, xpt)
      qfd = MATMUL_D(xpt, xptd, x, xd, qf)
    END SUBROUTINE QUADRATICFORM_D
    

    Test code is given below

    program main
    	implicit none
    	real :: A(2,2), AD(2,2), x(2), xd(2), qf, qfd
    
    	A = reshape( (/ 1,2,2,1 /), (/2,2/) )
    	AD = reshape( (/ -3,4,5,1 /), (/2,2/) )
    	x = (/ 3,4 /)
    	xd = (/ 7,2 /)
    
    	call quadraticForm(A,x,QF)
    
    	print *, A
    	print *, AD
    	print *, qf
    
    	! SUBROUTINE QUADRATICFORM_D(a, ad, x, xd, qf, qfd)
    	call QUADRATICFORM_D(A, AD, x, xd, qf, qfd)
    
    	print *, qfd
    
    end program main
    

    Test inputs are

    \mathbf{A}=\left[\begin{array}{cc}1 & 2\\2 & 1\end{array}\right]

    d\mathbf{A}=\left[\begin{array}{cc}-3 & 5\\4 & 1\end{array}\right]

    \mathbf{x}=\left[\begin{array}{c}3\\4\end{array}\right]

    d\mathbf{x}=\left[\begin{array}{c}7\\2\end{array}\right]

    and outputs are

     A  :    1.0000000       2.0000000       2.0000000       1.0000000
     AD :   -3.0000000       4.0000000       5.0000000       1.0000000
     QF :    73.000000
     dQF:    291.00000
    

    Manual calculation (not shown here) results in 291, also.

    Compilation and extra code is given below

    rem compilation of MATMUL example
    gfortran main.f90  quadraticForm.f90 quadraticform_d.f90 missingfunctions.f90   -o main.exe
    
    ! FILE: MissingFunctions.f90
    function matmul_d(a,ad,b,bd,ab)
    	implicit none
    	real, intent(in), dimension(:,:) :: a,b,ad,bd
    	real, intent(out), dimension( size(a,1), size(b,2) ) :: ab
    
    	real matmul_d( size(a,1), size(b,2) )
    
    	ab = matmul(a,b)
    	matmul_d = matmul(ad,b) + matmul(a,bd)
    end function matmul_d
    
    real function dot_product_d(a,ad,b,bd,ab)
    	implicit none
    	real, intent(in), dimension(:) ::  a,ad,b,bd
    	real, intent(out) :: ab
    	ab = dot_product(a,b)
    	dot_product_d = dot_product(ad,b) + dot_product(a,bd)
    end function dot_product_d
    
    function transpose_d(v, vd, vt)
    	implicit none
    	real, dimension(:,:), intent(in) :: v, vd
    	real, dimension(size(v,2), size(v,1) ) :: transpose_d
    	real, dimension(size(v,2), size(v,1) ), intent(out) :: vt
    
    	vt = transpose(v)
    	transpose_d = transpose(vd)
    end function transpose_d
    

    Notes:

    I have obtained strange results (automatically declared local variables of wrong size) from TAPENADE, when not used local variables for transposed variables. I don’t know the reason, but such behaviour of TAPENADE can avoided by some minor revisions to available code, as in the case of this example.

    An example for callbacks in ctypes

    November 15, 2009

    This example implements a rectangle (Riemann sum) quadrature rule in c compiled to a dll and used in Python via ctypes. File “quad.c” contains and exports

  • quadrature routine double rectquad(realfunction F,double x1,double x2,int N)
  • 2 sample functions in c to be integrated samplefun1 x \rightarrow x+2 x^2 and samplefun2 x \rightarrow x+n x^n
  • an integer variable npow as a parameter to one of the integrand function
  • Quadrature routine computes the following approximation

    \displaystyle \int_{x_{1}}^{x_{2}}f\left( x\right) dx\simeq \sum_{i=0}^{n}\left( \frac{x_{2}-x_{1}}{n}\right) f\left( x_{1}+i\frac{x_{2}-x_{1}}{n}\right)

    /* FILE quad.c */
    #include "dll.h"
    
    #ifdef __cplusplus
    extern "C" {
    #endif
    
    EXPORT int npow;
    
    typedef double realfunction(double x);
    
    EXPORT double rectquad(realfunction F,double x1,double x2,int N){
        double dA;
        double result = 0.0;
        double dx = (x2-x1)/N;
        double x = x1;
        int i;
        for(i=0;i<N;i++){
            dA = F(x+dx/2)*dx;
            result += dA;
            x += dx;
        }
        return result;
    }
    
    EXPORT double samplefun1(double x){
    	return x+2*x*x;
    }
    
    EXPORT double samplefun2(double x){
    	return x+npow*pow(x,npow);
    }
    
    #ifdef __cplusplus
    }
    #endif
    

    Code given above compiled to a dll by

    gcc -c quad.c -DBUILD_DLL
    gcc -shared -o quad.dll quad.o -Wl,--out-implib,libquad.a
    

    Generated dll file is used in following Python script. Quadrature routine is called with Python functions for x \rightarrow x and x \rightarrow x^2

    from ctypes import *
    import math
    
    quadlib = cdll.quad
    
    samplelibfun1 = quadlib.samplefun1
    samplelibfun2 = quadlib.samplefun2
    
    samplelibfun1.restype = c_double
    samplelibfun1.argtypes = [c_double]
    
    samplelibfun2.restype = c_double
    samplelibfun2.argtypes = [c_double]
    
    realFunction = CFUNCTYPE( c_double, c_double )
    
    npow = c_int.in_dll(quadlib, "npow")
    x = 2.0
    for i in range(10):
        npow.value = i
        print "npow:%i => result(c-dll):%12.7f, result(Python):%12.7f" % (npow.value, samplelibfun2( x ), x+i*math.pow(x,i) )
    
    quadlib.rectquad.argtypes = [realFunction, c_double, c_double, c_int]
    quadlib.rectquad.restype = c_double
    
    
    def Id(x):
        #print "x=",x
        return x
    
    def sqr(x):
        #print "x=",x
        return x*x
    
    print "_"*30
    print "quadrature of c function samplelibfun1"
    result = quadlib.rectquad( cast(samplelibfun1,realFunction ), c_double(0.), c_double(1.), c_int(64) )
    print "result=",result
    
    
    print "_"*30
    print "quadrature of c function samplelibfun2 with parameter npow=3"
    npow.value = 3
    result = quadlib.rectquad( cast(samplelibfun2,realFunction ), c_double(0.), c_double(1.), c_int(64) )
    print "result=",result
    
    
    print "_"*30
    print "quad of x->x"
    for n in range(25,150,25):
        result = quadlib.rectquad( realFunction(Id), c_double(0.), c_double(1.), c_int(n) )
        print "result for subdivision %3i is %12.8f" % (n, result)
    
    print "_"*30
    print "quad of x->x*x"
    for n in range(50,650,50):
        result = quadlib.rectquad( realFunction(sqr), c_double(0.), c_double(1.), c_int(n) )
        print "result for subdivision %3i is %12.8f" % (n, result)
    
    IDLE 2.6.2      ==== No Subprocess ====
    >>> 
    npow:0 => result(c-dll):   2.0000000, result(Python):   2.0000000
    npow:1 => result(c-dll):   4.0000000, result(Python):   4.0000000
    npow:2 => result(c-dll):  10.0000000, result(Python):  10.0000000
    npow:3 => result(c-dll):  26.0000000, result(Python):  26.0000000
    npow:4 => result(c-dll):  66.0000000, result(Python):  66.0000000
    npow:5 => result(c-dll): 162.0000000, result(Python): 162.0000000
    npow:6 => result(c-dll): 386.0000000, result(Python): 386.0000000
    npow:7 => result(c-dll): 898.0000000, result(Python): 898.0000000
    npow:8 => result(c-dll):2050.0000000, result(Python):2050.0000000
    npow:9 => result(c-dll):4610.0000000, result(Python):4610.0000000
    ______________________________
    quadrature of c function samplelibfun1
    result= 1.16662597656
    ______________________________
    quadrature of c function samplelibfun2 with parameter npow=3
    result= 1.24990844727
    ______________________________
    quad of x->x
    result for subdivision  25 is   0.50000000
    result for subdivision  50 is   0.50000000
    result for subdivision  75 is   0.50000000
    result for subdivision 100 is   0.50000000
    result for subdivision 125 is   0.50000000
    ______________________________
    quad of x->x*x
    result for subdivision  50 is   0.33330000
    result for subdivision 100 is   0.33332500
    result for subdivision 150 is   0.33332963
    result for subdivision 200 is   0.33333125
    result for subdivision 250 is   0.33333200
    result for subdivision 300 is   0.33333241
    result for subdivision 350 is   0.33333265
    result for subdivision 400 is   0.33333281
    result for subdivision 450 is   0.33333292
    result for subdivision 500 is   0.33333300
    result for subdivision 550 is   0.33333306
    result for subdivision 600 is   0.33333310
    >>> 
    

    Compiling a C code to a dll (MinGW) and using in Python via ctypes

    November 14, 2009

    In this example a dll file that computes the state derivative function

    \displaystyle \frac{dx_1(t)}{dt}=-(x_1^2+1) \\ \frac{dx_2(t)}{dt}=t \sin(x_1+x_2) \\ \frac{dx_3(t)}{dt}=\sin(\cos(x_1)+\sin(x_3))

    is generated and used in a Python script by ctypes module. Source files for the dll file are given below

    /* FILE dll.h */
    #ifdef BUILD_DLL
    /* DLL export */
    #define EXPORT __declspec(dllexport)
    #else
    /* EXE import */
    #define EXPORT __declspec(dllimport)
    #endif
    
    /* FILE statederivative.c */
    #include "dll.h"
    
    #ifdef __cplusplus
    extern "C" {
    #endif
    
    EXPORT void statederivative(double t,double x[], double dx[]){
        dx[0] = -(x[0]*x[0]+1);
        dx[1] = t*sin(x[0]+x[1]);
        dx[2] = sin( cos(x[0])+sin(x[2]) );
    }
    
    #ifdef __cplusplus
    }
    #endif
    

    If MinGW is used, dll file is created issuing following commands

    gcc -c statederivative.c -DBUILD_DLL
    gcc -shared -o statederivative.dll statederivative.o -Wl,--out-implib,libstatederivative.a
    
  • 1st line creates statederivative.o with define “BUILD_DLL” so that “__declspec(dllexport)” is tagged to the function declaration
  • 2nd line links the object file “statederivative.o” to generate “statederivative.dll” passing the link directives “–out-implib,libstatederivative.a” to linker to generate an import library.
  • To examine contents of the object code and dll use nm command

    >nm statederivative.o
    00000000 b .bss
    00000000 d .data
    00000000 i .drectve
    00000000 t .text
             U _sin
    00000000 T _statederivative
    
    > nm statederivative.dll
    67c44074 b .bss
    67c41170 t .text
    ...
    67c44010 b _next_atexit
    67c415f8 T _sin
    67c41180 T _statederivative
    67c43000 r _w32_atom_suffix
    ...
    

    Following script illustrates how to use the generated dll file “statederivative.dll”

    from ctypes import *
    import math
    
    sdlib = cdll.statederivative
    print "sdlib:",sdlib
    twodoubles = c_double*2
    threedoubles = c_double*3
    
    print "_"*30
    print "twodoubles  :",twodoubles
    print "threedoubles:",threedoubles
    
    print "_"*30
    t = c_double(1.0)
    x = threedoubles(3., 4., 5.)
    dx = threedoubles()
    dxPython = threedoubles()
    
    print "dx[0],dx[1],dx[2]:",dx[0],dx[1],dx[2]
    print " x[0], x[1], x[2]:",x[0],x[1],x[2]
    
    # declare signature for function in dll
    sdlib.statederivative.argtypes = [c_double, threedoubles, threedoubles]
    
    # define C-dll function in Python
    def statederivative(t,x,dx):
        dx[0] = -(x[0]*x[0]+1)
        dx[1] = t*math.sin(x[0]+x[1])
        dx[2] = math.sin( math.cos(x[0])+math.sin(x[2]) )
    
    # call state derivative routine in dll
    sdlib.statederivative( t, x, dx)
    # call state derivative Python implementation
    statederivative(t.value, x, dxPython)
    
    # print dll and Python results
    print "_"*30
    print "### DLL-COMPUTATION ###"
    print " x[0], x[1], x[2]:",x[0],x[1],x[2]
    print "dx[0],dx[1],dx[2]:",dx[0],dx[1],dx[2]
    print "### PYTHON-COMPUTATION ###"
    print "dx[0],dx[1],dx[2]:",dxPython[0],dxPython[1],dxPython[2]
    
    # perform Euler integration using c-dll function
    print "_"*30
    print "Euler integration"
    h = 0.25
    t = 0
    x[0] = 0
    x[1] = 2
    x[2] = 10
    solution = [[t,x[0],x[1],x[2]]]
    for i in range(30):
        sdlib.statederivative( c_double(t), x, dx )
        x[0] += ( h * dx[0] )
        x[1] += ( h * dx[1] )
        x[2] += ( h * dx[2] )
        t += h
        print "(t:%10.6f)\tx[0:2]:%10.6f,%10.6f,%10.6f" % (t,x[0],x[1],x[2])
    
        solution.append( [t,x[0],x[1],x[2]] )
    
    solution = sp.array( solution ) # convert to list of lists to matrix
    print solution
    
    
  • line 7 declares a ctypes array of double declaration “threedoubles = c_double*3
  • line 23 declares c function signature using ctypes built-in c types and the type “threedoubles" declared in script
  • Output of this Python script in given below

    IDLE 2.6.2      ==== No Subprocess ====
    >>>
    sdlib: <CDLL 'statederivative', handle 67c40000 at baa150>
    ______________________________
    twodoubles  : <class '__main__.c_double_Array_2'>
    threedoubles: <class '__main__.c_double_Array_3'>
    ______________________________
    dx[0],dx[1],dx[2]: 0.0 0.0 0.0
     x[0], x[1], x[2]: 3.0 4.0 5.0
    ______________________________
    ### DLL-COMPUTATION ###
     x[0], x[1], x[2]: 3.0 4.0 5.0
    dx[0],dx[1],dx[2]: -28.0 0.656986598719 -0.929360160426
    ### PYTHON-COMPUTATION ###
    dx[0],dx[1],dx[2]: -28.0 0.656986598719 -0.929360160426
    ______________________________
    Euler integration
    (t:  0.250000)	x[0:2]: -0.250000,  2.000000, 10.110085
    (t:  0.500000)	x[0:2]: -0.496094,  2.061499, 10.192514
    (t:  0.750000)	x[0:2]: -0.715570,  2.186497, 10.238486
    ...
    (t:  7.000000)	x[0:2]: -1.000000,  4.141593,  9.996768
    (t:  7.250000)	x[0:2]: -1.000000,  4.141593,  9.996517
    (t:  7.500000)	x[0:2]: -1.000000,  4.141593,  9.996319
    [[  0.           0.           2.          10.        ]
     [  0.25        -0.25         2.          10.11008535]
     [  0.5         -0.49609375   2.06149912  10.19251403]
    ...
    
     [  7.          -1.           4.14159265   9.99676793]
     [  7.25        -1.           4.14159265   9.99651692]
     [  7.5         -1.           4.14159265   9.99631868]]
    >>>
    

    lines 13 and 15 give the c and Python results of statederivative function.
    In following Python script the state derivative function integrated with one of SciPy ODE-IVP solvers/integrators. For this purpose, a Python wrapper function “sd_scipy_wrapper(y,t)” is implemented and given as input argument in SciPy function “spint.odeint(y,t)“. Finally SciPy integration and previously obtained Euler integration results are plotted using matplotlib. Note that there is some overhead due to using wrapper routine

    import scipy as sp
    import scipy.integrate as spint
    import pylab as pl
    
    solution = sp.array( solution ) # convert to list of lists to matrix
    
    def sd_scipy_wrapper(y,t):
        dy = threedoubles()
        sdlib.statederivative(c_double(t),threedoubles(*y),dy)
        result = [dy[i] for i in range(3)]
        print result
        return [dy[i] for i in range(3)]
    
    T = sp.arange(0,10,0.1)
    Y0 = [0.,2.,10.]
    Y = spint.odeint(sd_scipy_wrapper, Y0, T)
    
    fig1 = pl.figure()
    fs = 12
    pl.subplot(311)
    pl.plot(Y[:,0],Y[:,1])
    pl.xlabel("$X_1$", size=fs)
    pl.ylabel("$X_2$", size=fs)
    pl.grid()
    
    pl.subplot(312)
    pl.plot(Y[:,0],Y[:,2])    pl.xlabel("$X_1$", size=fs)
    pl.ylabel("$X_3$", size=fs)
    pl.grid()
    
    pl.subplot(313)
    pl.plot(Y[:,1],Y[:,2])
    pl.xlabel("$X_2$", size=fs)
    pl.ylabel("$X_3$", size=fs)
    pl.grid()
    
    fig2 = pl.figure()
    pl.plot(T,Y[:,0],T,Y[:,1],T,Y[:,2])
    fig2.hold('on')
    pl.plot(solution[:,0],solution[:,1],'o', solution[:,0],solution[:,2],'o', solution[:,0],solution[:,3], marker='o',markeredgecolor='None')
    leg = pl.legend(["$X_1 [\mathtt{SciPy}]$","$X_2 [\mathtt{SciPy}]$","$X_3 [\mathtt{SciPy}]$","$X_1 [\mathsf{Euler}]$","$X_2 [\mathsf{Euler}]$","$X_3 [\mathsf{Euler}]$"])
    leg.get_frame().set_alpha(0.5)
    pl.xlabel("$T$")
    pl.ylabel("$X_j$")
    pl.grid()
    
    pl.show()
    
    fig1.savefig('fig1.png')
    fig2.savefig('fig2.png')
    

    Solution time histories for example systemState-space phase portraits

    Basic F2PY Examples (windows, MinGW)

    November 12, 2009

    To extend Python with following FORTRAN subroutine that computes the \left(x,y\right)\rightarrow \sin(x+y)

    subroutine f1(x,y,z)
    	real, intent(in) :: x,y
    	real, intent(out) :: z
    	
    	z = sin(x+y)
    end subroutine f1
    

    f2py is executed as follows

    rem call f2py
    f2py.py --build-dir .\tmp --fcompiler=gnu95 --compiler=mingw32 -c lib1.f90 -m lib1
    

    Command line arguments above indicate that

  • temporary directory for generated code is .\tmp
  • fortran compiler is g95 and c compiler/linker is mingw32
  • build extension modules from source file(s) lib1.f90
  • name of generated library is lib1
  • To determine supported compilers, type

    f2py.py -c --help-fcompiler
    

    In my case F2PY found 2 compilers

    ...
    Fortran compilers found:
      --fcompiler=gnu    GNU Fortran 77 compiler (3.4.5)
      --fcompiler=gnu95  GNU Fortran 95 compiler (4.5.0)
    ....
    

    Upon completion of F2PY call, lib1.pyd is generated. Example usage script for generated module is given next

    import lib1
    
    print "_"*30
    print "Contents of lib1"
    print dir(lib1)
    print lib1.__doc__
    
    print "_"*30
    print "Help for lib1 function f1"
    print lib1.f1.__doc__
    
    print "_"*30
    print "calling F2PY generated extension module function"
    z = lib1.f1(1., 2.)
    print "lib1.f1(1., 2.) = ", z
    import math
    print "     sin(1.+2.) = ", math.sin(1.+2.)
    
    IDLE 2.6.2      ==== No Subprocess ====
    >>> 
    ______________________________
    Contents of lib1
    ['__doc__', '__file__', '__name__', '__package__', '__version__', 'f1']
    This module 'lib1' is auto-generated with f2py (version:2_6832).
    Functions:
      z = f1(x,y)
    .
    ______________________________
    Help for lib1 function f1
    f1 - Function signature:
      z = f1(x,y)
    Required arguments:
      x : input float
      y : input float
    Return objects:
      z : float
    
    ______________________________
    calling F2PY generated extension module function
    lib1.f1(1., 2.) =  0.141120001674
         sin(1.+2.) =  0.14112000806
    >>>