home *** CD-ROM | disk | FTP | other *** search
/ Hall of Fame / HallofFameCDROM.cdr / oilfield / dyngraf2.lzh / SPEDYN.BAS < prev    next >
BASIC Source File  |  1989-03-09  |  37KB  |  1,067 lines
  1. 'University of Oklahoma Sucker Rod Pumping Unit Diagnostic System
  2. '
  3. 'By Robert P. Moore
  4. '
  5. 'The Diagnostic system is designed to evaluate surface dynagrphs
  6. 'and generate the pump dynagraph from the input data.
  7. '
  8. 'The diagnostic system was developed in 1988 - 89 as partial
  9. 'fulfilment of the requirements for a Master's of Science in
  10. 'Petroleum Engineering.  The system is designed for operation on an
  11. 'IBM Compatable PC which is attached to a line printer, Calcomp
  12. 'Model 23120 12 button Digitizer, and a pen plotter.
  13.  
  14. 'Both the digitizer and the pen plotter are configured through
  15. 'COM1.  An A-B switch is needed to select which component is needed
  16. 'at a particular time.  The printer is attached to the parallel
  17. 'interface.  An Intel 8087-2 Math Co-Processor was also installed to
  18. 'speed operations.
  19.  
  20. 'The programming language is Borland's Turbo BASIC for the program
  21. 'and Hewlett Packard Graphics Language for the plotter.  Also
  22. 'included on the disc is a LOTUS - 123 spreadsheet which is used to
  23. 'plot the data if a pen plotter is not available.  The program
  24. 'creates a sequential data file which can be loaded into the spread-
  25. 'sheet.  A macro creates and saves the graphs for plotting with the
  26. 'PRINTGRAPH utility.  A list of input and output is given below.
  27.  
  28. 'The purpose of the dynagraph diagnostic package is to enable a
  29. 'production engineer to analyse a dynagraph quickly and easily
  30. 'through a PC environment.  Total run time is less than five minutes
  31. 'in all cases evaluated to date, including operator input.
  32. 'The future in this area could be portable PC's and dynagraph
  33. 'interfaces which allow dynagraph analysis on-site and almost in
  34. 'real time.
  35.  
  36. 'A field dynagraph is included with this disk to enable operation of
  37. 'the program.  Also, output generated is given.
  38.  
  39. 'The author would like to recognize the previous work of Cox and
  40. 'Chacin.  A modified version of Cox's dynagraph input routine is
  41. 'while the downhole dynagraph generation solution and routine were
  42. 'developed by Chacin as part of his Doctoral research.  I would also
  43. 'like to thank Dr. Ron Evans for permission to submit this program.
  44. '
  45. '
  46. '.....  Input Data....
  47.  
  48. ' 1. Well Name
  49. ' 2. Surface Dynagraph Profile
  50. ' 3. Peak Polished Rod Load
  51. ' 4. Minimum Polished Rod Load
  52. ' 5. Travelling Valve Load
  53. ' 6. Standing Valve Load
  54. ' 7. Counter Balance Effect
  55. ' 8. Rod Strin Data
  56. ' 9. Fluid Gravity (Optional)
  57. '10. Damping Factor
  58. '11. # Strokes / Minute
  59. '12. Stroke Length
  60. '13. Pump Plunger Diameter
  61. '14. Tubing Head Pressure
  62.  
  63.  
  64. '.... Ouput Data ....
  65.  
  66. ' 1. Calculated Fluid Gravity
  67. ' 2. Formation Flowing Pressure
  68. ' 3. Workin Fluid Depth
  69. ' 4. Maximum Polished Rod Load
  70. ' 5. Maximum Rod Stress in Top Rod Section
  71. ' 6. Average Maximum Polished Rod Load
  72. ' 7. Average Minimum Polished Rod Load
  73. ' 8. Ideal Counter Balance Effect
  74. ' 9. Polished Rod Horsepower
  75. '10. Maximum Pumpstroke
  76. '11. Effective Pumpstroke
  77. '12. Theoretical Pumping Rate
  78. '13. Peak Torque
  79. '14. Peak Torque Crank Angle
  80. '15. Plot of Pump Load vs. Pump Position
  81. '16. Plot of Torque vs. Crank Angle
  82.  
  83. 'If you have any questions I can be reached at 364-8137.
  84.  
  85. '... Digitize the Dynagraph ...
  86.  
  87. $com1 4096                'Dimension Com Port Buffer
  88. screen 2
  89. '
  90. begin:
  91. clear
  92. cls
  93.  
  94. line (0,0)-(0,0),3                      'Set up plot capability
  95. dim cr(20)                               '
  96. get(0,0)-(0,0),cr
  97. dim x(5000),y(5000),c$(5000)        'Dimension Arrays
  98. dim xg(5000),yg(5000)
  99.  
  100. input "Enter Well Name: ",wellid$
  101.  
  102. cls
  103. close #1
  104. open "com1:9600,n,8,1" as #1:out &H3FB,11     'Initialize tablet
  105. print #1, "2A"
  106. close #1
  107.     open "com1:9600,e,7,1"as #1          'Set tablet for point mode
  108.         print #1,"SQj2l"
  109. '
  110.     print "Digitize topmost point of dynagraph"
  111.             print #1,"?"
  112.             input #1,xymax,ymax,c$
  113.             print "Point digitized"
  114.             print "
  115.  
  116.     print "Digitize bottom-most point of dynagraph"
  117.             print #1,"?"
  118.             input #1,xymin,ymin,c$
  119.             print "Point digitized"
  120.     delay 1
  121.  
  122.         cls
  123.         print "Move cursor to beginning point of dynagraph and press "
  124.         print "button #0 when ready.  Digitize the upstroke by slowly"
  125.         print "tracing the upper half of the profile of the card."
  126.         print "
  127.         print "After digitizing the upstroke of the card,
  128.         print "press button #1 and move the cursor minutely to the right.
  129.         print "Release the button and trace the downstroke profile.
  130. '
  131. 430     print #1,"?"
  132.             input #1,xleft,yleft,c$
  133.                 if c$<>"1" then goto 430
  134. '
  135.         print #1,"SRj2l"        'Put tablet into Run-Prompt mode
  136.         i=0
  137.         print "           x       y       "
  138. '
  139. inputroutine:
  140. 410    i=i+1
  141.     j=i-1
  142. 420    print #1,"?"
  143.         input #1,x(i),y(i),c$(i)
  144.         if(x(i)<=x(j)) then goto 420
  145.             print "i=" i, x(i),y(i),c$(i)
  146.             if c$(i)="2" then goto 450
  147.         goto 410
  148. '
  149. 450     xright=x(i):yright=y(i):num=i
  150.  
  151. 455     i=i+1
  152.     j=i-1
  153. 460    print #1,"?"
  154.         input #1,x(i),y(i),c$(i)
  155.         if(x(i)>=x(j)) then goto 460
  156.             if (x(i)<=xleft) then goto 470
  157.             print "i=" i, x(i),y(i),c$(i)
  158.         goto 455
  159. '
  160. 470     locate 20,55
  161.         print "Digitizing is complete.
  162.         cnt = i
  163.         delay 2
  164. '
  165.  PlotCard:
  166.        cls
  167.        locate 1,30
  168.        print wellid$
  169.       xl=cint((xleft-2000)*.053)
  170.       yl=cint((yleft-6000)*-(.05/2))
  171. 'print xl,yl
  172. put(xl,yl),cr
  173.         for k=1 to i
  174.             xg(k)=cint((x(k)-2000)*.053)
  175.         yg(k)=cint((y(k)-6000)*-(.05/2))
  176.             line -(xg(k),yg(k))
  177.         next k
  178. '
  179.         locate 20,1
  180.         input "Was the digitized shape correct?(y/n)";yes$
  181.         if (yes$="Y") then goto start
  182.         if (yes$="y") then goto start
  183.         print
  184.         input "Redo the digitizing. Press enter to start again.";yes$
  185.         goto begin
  186.  
  187. start:
  188. cls
  189. '
  190. input "Enter value for Peak Polished Rod Load";lodmx
  191. print "
  192. input "Enter value for Minimum Polished Rod Load";loadmin
  193. print "
  194. input "Enter Travelling Valve Test Value";tv%
  195. print "
  196. input "Enter Standing Valve Test Value";sv%
  197. print "
  198. input "Enter Counter Balance Effect Value";cbe%
  199. print "
  200.  
  201.                         'Dimension Arrays
  202. dim posi(5000),load(5000)
  203. dim posit(5000),newld(5000)
  204. '
  205. i = cnt
  206.  
  207. x(i)=x(1)                   'Convert data to dynagraph scale    
  208. y(i)=y(1)
  209. '
  210. for j=1 to i
  211.        posi(j)=ceil((x(j)-xleft)*(1000/(xright-xleft)))
  212.            load(j)=ceil(((y(j)-ymax)*((lodmx-loadmin)/(ymax-ymin)))+lodmx)
  213.            if(posi(j-1)=1000) and (posi(j)=1000) then goto 508
  214. '          print posi(j),load(j)
  215. 508   next j
  216. '
  217. print ""
  218. print "Data Values Converted to Desired Load and Position Values......Wait"
  219. print "
  220. '
  221. EliminateDuplicates:
  222.     posi(1)=0
  223.        n=0                                'Eliminate duplicate x-axis points
  224. 590     n=n+1
  225.         m=n+1
  226.         if (posi(n)>=1000) then posi(n)=1000: goto 591  'Eliminate duplicates
  227.         if (posi(m)<=posi(n)) then posi(m)=(posi(n))+1   'on upstroke
  228.         goto 590
  229. '                                                       'Eliminate duplicates
  230. 591     if (posi(m)>=posi(n)) then posi(m)=posi(n)-1    'on downstroke
  231.     if (posi(n)=0) then print "Duplicate Points Eliminated..Wait":goto interpolate
  232.         n=n+1
  233.         m=n+1
  234.         goto 591
  235. '
  236. interpolate:                    'Interpolate between data
  237.     numb=n                                  'to obtain 1000 data points
  238.     posmax=posi(numb)                       'on the upstroke and the
  239.     n=0:m=0:j=0                             'downstroke.
  240. 502    j=1
  241.     n=1
  242. 505     m=n+1
  243.     if posi(n)=1000 then goto 592
  244.         pt1=posi(n):ld1=load(n)
  245.     pt2=posi(m):ld2=load(m)
  246. 510    if((pt2-pt1)=1) then goto 520
  247. 515    newpt2=pt1+ceil((pt2-pt1)/2)
  248.     newload=ld1+ceil((ld2-ld1)/2):goto 535
  249. '
  250. 520     posit(j)=pt1
  251.     newld(j)=ld1
  252.     a=1
  253.     goto 530
  254. '
  255. 535    posit(j)=pt1:newld(j)=ld1
  256.     posit(j+1)=newpt2:newld(j+1)=newload
  257.     a=2
  258. 530    j=j+a
  259.         n=n+1
  260.         goto 505
  261. '
  262. 525     m=n+1
  263. 592     pt1=posi(n):ld1=load(n)
  264.     pt2=posi(m):ld2=load(m)
  265. 593    if((pt1-pt2)=1) then goto 595
  266. 594    newpt2=pt1-ceil((pt1-pt2)/2)
  267.     newload=ld1-ceil((ld1-ld2)/2)
  268.     goto 597
  269. '
  270. 595     posit(j)=pt1
  271.     newld(j)=ld1
  272.     a=1
  273.     goto 596
  274. '
  275. 597    posit(j)=pt1:newld(j)=ld1
  276.     posit(j+1)=newpt2:newld(j+1)=newload
  277.     a=2
  278. 596     if (posit(j)=0) then goto 540
  279.     j=j+a
  280.         n=n+1
  281.         goto 525
  282. '
  283. 540    print "num=" j
  284.     if (j>2001) then print "Error in Interpolation":end
  285. repeat:
  286.     for b=1 to j
  287.         posi(b)=posit(b)
  288.         load(b)=newld(b)
  289.     next b
  290.     if(j<2001) then goto 502
  291. '
  292.     print "Interpolation Completed"
  293.     print "
  294. '
  295.  
  296. Filter:
  297.         b = 0
  298.         for a% = 1 to 500
  299.              y = cint((2001/501)*a%)
  300.              posi(a%) = posi(y)
  301.              load(a%) = load(y)
  302. '            print a%,posi(a%),load(a%)
  303.         next a%
  304. 'lprint "Filter",a%,posi(a%),load(a%)
  305. cls
  306.  
  307.  
  308.    lprint "     *********************************"
  309.    lprint
  310.    lprint "     "Wellid$:lprint
  311.    lprint "     *********************************"
  312.    lprint
  313.  
  314.    input "Enter number of different rod sections:",ntap%
  315.    lprint "     Rod String Data":lprint:lprint
  316.  
  317.      '...Dimension and initialize variables ...
  318.  
  319.                      'ntap% = # of tapers
  320.    Dim irod%(ntap%)  'rod material ; 1 = steel; 2 = fiberglass
  321.    dim dia(ntap%)    'diameter in inches of the rods
  322.    dim aleng(ntap%)  'length of rod section,input in ft converted to in
  323.    dim aaleng(ntap%) 'part of aaleng that will need to be interpolated
  324.    dim slfact(ntap%) 'strain multiplier to give load
  325.    dim wr(ntap%+1)   'weight of rods
  326.    dim ey(ntap%)     'young's modulus
  327.    dim a(ntap%)      'speed of wave propogation in the rod section
  328.    dim amc(ntap%)    'weight per foot of the rod section -- API RP-11l
  329.    totdepth = 0
  330.  
  331.      '... Input data from screen ...
  332.    print:print"Enter rod data top to bottom.":print
  333.    for i% = 1 to ntap%
  334.    print:print"Rod section number",i%
  335.    print:input"Enter diameter of this rod section #.### in :",dia(i%)
  336.    print:input"enter length of this rod section - ft :",aleng(i%)
  337.          totdepth = totdepth + aleng (i%)
  338. '         print totdepth
  339.          aleng(i%) = aleng(i%) * 12      'convert to inches
  340. 1  print:input"Enter rod type (1 = steel, 2 = fiberglass) :",irod%(i%)
  341.    if irod%(i%) < 1 or irod%(i%) > 2 then print "bad input" : goto 1
  342.    lprint "       Rod Section Number "i%:
  343.    if irod%(i%) = 1 then rodtype$ = "steel"
  344.    if irod%(i%) = 2 then rodtype$ = "fiberglass"
  345.    lprint "         "aleng(i%)/12 " feet of "dia(i%) " inch " rodtype$ " rods."
  346.    lprint
  347.    next i%
  348.  
  349.         '... Input fluid gravity or calculate it from the SV check ...
  350.  
  351.    lprint "     *********************************"
  352.    lprint
  353.  
  354.    res = 5!         'convergence criterion on the position variable
  355.    cut = 1!         'maximum length ignored for interpolation
  356.    cut = cut * 12
  357.    pi = 3.1415926#
  358.  
  359.      '... dimension variables ...
  360.  
  361.    m% = 500         '500 data points will be read
  362.    dim time(m%+1)
  363.    dim ww(m%+1), uprime(ntap%+1,m%+1),vprime(ntap%+1,m%+1),wprime(ntap%+1,m%+1)
  364.    dim us(2*m%+1),vs(2*m%+1),ws(2*m%+1)
  365.                     'in the previously dimensioned variables, u refers to
  366.                     'position in inches, v refers to strain, and w refers
  367.                     'to velocity in in / sec, time is in seconds
  368.  
  369.      ' ... Initialize arrays ...
  370.  
  371.   ' for i% = 1 to m%
  372.   '      time(i%+1) = time(i%) + adt
  373.   ' next i%
  374.  
  375.    for i% = 1 to ntap% + 1
  376.         wr(i%) = 0
  377.    next i%
  378.  
  379.      ' ...Loop to calculate rod section properties ...
  380.  
  381.    for i% = 1 to ntap%
  382.       if irod%(i%) = 2 then goto 10
  383.           ' ... for steel rods ...
  384.           a(i%)  = 195600!
  385.           ey(i%) = 3.05E+07
  386.  
  387.           if dia(i%) = 1.5   then amc(i%) = 533!*pi*(dia(i%)^2)/6912!  'lb/in
  388.           if dia(i%) = 1.25  then amc(i%) = 4.538/12!
  389.           if dia(i%) = 1.125 then amc(i%) = 3.676/12!
  390.           if dia(i%) = 1!    then amc(i%) = 2.904/12!
  391.           if dia(i%) = .875  then amc(i%) = 2.224/12!
  392.           if dia(i%) = .75   then amc(i%) = 1.634/12!
  393.           if dia(i%) = .625  then amc(i%) = 1.135/12!
  394.           if dia(i%) = .5    then amc(i%) = 0.726/12!
  395.           goto 20
  396.  
  397. 10     '   ... for fiberglass rods ...
  398.        a(i%)       = 177270!
  399.        ey(i%)      = 7200000!
  400.        amc(i%)     = 153! * pi * (dia(i%)^2)/(4! * 144! * 12!)
  401.  
  402. 20     slfact(i%)  = pi * (dia(i%)^2) * ey(i%) / 4!
  403.        wr(i%)      = amc(i%) * aleng(i%)
  404.  
  405.        'print "20","i%="i%,"dia(i%)="dia(i%),"aleng(i%)=",aleng(i%)
  406.        'print "a(i%)="a(i%),"ey(i%)="ey(i%),"slfact(i%)="slfact(i%)
  407.        'print "wr(i%)="wr(i%)
  408.    next i%
  409.  
  410.       ' ... calculate total rod weight and volume ...
  411.  
  412.    twr    = 0!
  413.    rodvol = 0!
  414.  
  415.    for i% = 1 to ntap%
  416.         rodvol = rodvol + ( pi * (dia(i%)^2)/4!) * aleng(i%)
  417.         twr    = twr + wr(i%)
  418.    next i%
  419.  
  420.  
  421.    '   ... input dynagraph data ...
  422.  
  423.  
  424.       lprint "     Recorded Travelling Valve Load = "tv%" pounds":lprint
  425.       lprint "     Recorded Standing Valve Load = "sv%" pounds":lprint
  426.       lprint "     Recorded Counter Balance Effect = "cbe%" pounds":lprint
  427.       lprint "     *********************************":lprint:lprint
  428.  
  429.      if sv% = 0 then goto 5   'sv load is not input, fluid gravity must be
  430.      print
  431.      input "Do you want to input a fluid gravity instead of allowing a value be calculated?  ",grvans$
  432.        if grvans$ = "y" then goto 5
  433.        if grvans$ = "Y" then goto 5 else goto 7
  434.  
  435. 5  print:input  "Specific gravity of fluid :",grav
  436.    lprint using "     Input Fluid Gravity =#.###";grav:lprint
  437.    goto 8
  438.  
  439. 7       grav = 7.874 * ( 1 - sv% / twr )        ' calculate fluid gravity
  440.         lprint using"     Calculated Fluid Gravity =#.#### "; grav:lprint
  441.  
  442. 8  print:input  "Enter damping factor (0.1 1/sec is recommended) :",c
  443.    lprint using"     Damping factor (1/sec) =#.###";c:lprint
  444.  
  445.  
  446.  
  447.       lprint
  448.       lprint "     *********************************"
  449.       lprint:lprint
  450. ' ... calculate buoyed weight of rods - wrf ...
  451.  
  452.    buo = rodvol * grav * 62.4 / 1728!
  453.    wrf = twr - buo
  454.  
  455.      ' ... Calculate rod stretch ...
  456.  
  457.    rodst = 0
  458.    for i% = 1 to ntap%
  459.         rodst = rodst + 32.2 * ((aleng(i%)/a(i%))^2)*(12!)
  460.         rodst = rodst + (4!*aleng(i%)/(pi*(dia(i%)^2)*ey(i%)))*(-buo+wr(i%+1))
  461.    next i%
  462.  
  463.    if ntap% = 3 then goto 30 else goto 40
  464.  
  465. 30 rodst = rodst + (4! * aleng(1) / (pi* (dia(1)^2)* ey(1))) * wr(3)
  466.  
  467. 40 if ntap% = 4 then goto 50 else goto 60
  468.  
  469. 50 rodst = rodst + (4! * aleng(1) / (pi* (dia(1)^2)* ey(1))) * (wr(3)+wr(4))
  470.    rodst = rodst + (4! * aleng(2) / (pi* (dia(2)^2)* ey(2))) * (wr(4))
  471.  
  472. 60 '
  473.  
  474. 'lprint "60","rodvol="rodvol,"twr="twr,"buo="buo,"wrf="wrf,"rodst="rodst
  475.  
  476.    ' ... Input unit data ...
  477.    lprint"     Pump Data":lprint
  478.    print: input "Input number of strokes / minute:",nst#
  479.    lprint "          Number of Strokes / Minute "nst#:lprint
  480.  
  481.    print: input "Input stroke Length - in:",stl
  482.    lprint "          Stroke Length= "stl" inches":lprint
  483.  
  484.    print: input "Input Pump Plunger Diameter - inches:",plsize
  485.    lprint "          Pump Plunger Diameter ="plsize"inches":lprint
  486.    print: input "Tubing Head Pressure - psi:",thp
  487.    lprint "          Tubing Head Pressure ="thp"psi":lprint
  488.  
  489.  
  490.         plunarea = pi / 4 *(plsize)^2
  491.  
  492.         if (sv% > 0 and tv% > 0) then goto 23 else goto 24
  493.  
  494. 23        pwf = int((twr - sv%) / (pi * (dia(ntap%)/2)^2) - (tv% - sv%) / plunarea + (thp * ( 1 - (pi * (dia(ntap%)/2)^2) / plunarea)))
  495.  
  496.    lprint "          Flowing Pressure ="pwf"psi":lprint
  497.           wflddp =  int(totdepth - (2.31 / grav) * ( pwf - thp))
  498.  
  499.    lprint "          Working Fluid Depth ="wflddp"feet":lprint
  500.  
  501. 24 lprint "     *********************************"
  502.    lprint:lprint
  503.  
  504.        's = stroke indicator s = 0 =>upstroke; s= 1 =>downstroke
  505.        'tst = timestep , changes by adt at every ct point
  506.        'ct  = counter for ouput of equally spaced (adt) position and strain
  507.        '      values
  508.  
  509.    s = 0 : tst = 0 : ct% = 1
  510.    upstarea = 0: downarea = 0
  511.    adt  = (60 / nst# / 500)                '500 dt points
  512.    'lprint "adt="adt
  513.    halfstr = stl /2                        '1/2 strokelength
  514.    dimlen% = 500                           '# of inputs from dynagraph file
  515.  
  516.       'tm(j) = time that the data was recorded,start of upstroke = 0
  517.       'posi(j) = position data from dynagraph file, converted to inches
  518.       '          varying from - stl/2 to + stl/2
  519.       'load(j) = load value from dynagraph file, converted to strain
  520.  
  521.    '   ... dimension, load, and convert dynagraph data to time basis ...
  522.  
  523.    '            This program is set up to use a dynagraph file with
  524.    '            250 points on the upstroke and downstroke.  If a
  525.    '            different size file is used the program will need to be changed.
  526.  
  527.    dim tm(dimlen%),sload(dimlen%)
  528.    maxload = 0                    'Maximum load in first rod section
  529.  
  530.    '... Read load and position data from file ...
  531.  
  532.         for j% = 1 to 500
  533.  
  534.              posi(j%)=(posi(j%)/ 1000 * stl)-halfstr 'sine wave w/ amp = stl/2
  535. '             lprint j%,posi(j%)
  536.              if j% > 250 then s = 1                 'downstroke
  537.  
  538.             '... Calculate area under the dynagraph using Trapezoidal Rule ...
  539.  
  540.                   if j% > 1 and j% < 500 then goto 61 else goto 62
  541. 61                     avload = load(j%) - ((load(j%) - load(j%-1))/2)
  542.                          if load(j%) > maxload then maxload = load(j%)
  543.                        area = avload * (posi(j%) - posi(j%-1))
  544.                        if s = 0 then upstarea = upstarea + area
  545.                        if s = 1 then downarea = downarea - area
  546.  
  547. 62           sload(j%) = (load(j%) - wrf)/slfact(1)   'load to strain
  548.              'print "pos=",posi(j%),"load=",sload(j%),"s=",s:delay .2
  549.  
  550.  
  551.                  sine   = posi(j%)/halfstr
  552.                  cosine = ((1-(sine)^2)^.5)
  553.  '                lprint sine,cosine
  554.                  if (cosine = 0 and posi(j%) < 0) then angle = -90  'avoid
  555.                  if (cosine = 0 and posi(j%) > 0) then angle =  90  'error
  556.  
  557.                  if cosine <> 0 then angle = (atn((sine)/(cosine)))*180 /pi
  558.                  'print sine,cosine,angle,posi(j%):delay 3
  559.                  tm(j%) =s*60/nst#+(-1)^s*((1/nst#)*60/4 + angle/(nst#*360)*60)
  560.                  'if tm(j%) < 0.5 then lprint j%,"time=",tm(j%),"posi=",posi(j%)
  561.         next j%
  562.  
  563.         lprint "      Dynagraph Data"
  564.  
  565.     maxprl  =  (upstarea) / stl
  566.     maxstres = maxload / (pi * ((dia(1))^2)/4)
  567.     lprint
  568.     lprint using"          Maximum Polished Rod Load (lb) =##### ";maxload:lprint
  569.     lprint using"          Max Rod Stress in Top String(psi) = ##### ";maxstres:lprint
  570.     lprint using"          Average Max Polished Rod Load (lb) =##### "; maxprl:lprint
  571.  
  572.     minprl  =  downarea / stl
  573.     lprint using"          Average Min Polished Rod Load (lb) =##### ";minprl
  574.     lprint
  575.  
  576.     idcbeff = .5 * ( maxprl + minprl )
  577.     lprint using"          Ideal Counter Balance Effect (lb)=##### "; idcbeff
  578.     lprint
  579.  
  580.     prhp    = stl * nst# * (upstarea - downarea) / ( 396000 * stl )
  581.      lprint using"          Polished Rod Horse Power =###.# "; prhp :lprint
  582.      lprint "     *********************************"
  583.  
  584.      lprint:lprint
  585.  
  586.  
  587.    '   ... Obtain 500 equally spaced(adt) data points ...
  588.    '             The strain and position values are interpolated to correspond
  589.    '             to the time step points
  590.    tm(1) = 0 : posi(1) = - stl/2
  591.    tm(500) = 60 / nst# : posi(500) = - stl/2
  592.    'chdir"c:\123\files"
  593.    'file$ = "erase2.prn"
  594.    'open file$ for output as #2
  595.    for k% = 1 to 500
  596.  
  597. 55           q = abs(tm(k%) - tst):'lprint ct%,"tst="tst,"tm="tm(k%),"decimal q = ",q
  598.              if ( q < .000001)then goto 56 else goto 58
  599. 56           time(ct%) = tst:'print time(ct%)
  600.              uprime(0,ct%) = - posi(k%):'lprint uprime(0,ct%)
  601.              vprime(0,ct%) = sload(k%):'lprint vprime(0,ct%)
  602.              goto 71
  603. 58           q = int(int(q*100) + 1)     '# of adt points between the values
  604.              'lprint ct%,"q=",q
  605.              if (tm(k%-1) < tst and tm(k%) >tst) then goto 59 else goto 65
  606.  
  607. 59           for i% = 1 to q
  608.              'lprint"i%",i%
  609.              time(ct%) = tst                            'interpolate strain and
  610.              intrp = - (tm(k%) - tst)/(tm(k%)-tm(k%-1)) 'position to time intv
  611.              uprime(0,ct%) = -(intrp*(posi(k%) - posi(k%-1)) + posi(k%))
  612.              vprime(0,ct%) = intrp*(sload(k%) - sload(k%-1)) + sload(k%)
  613. 70      'write #2,ct%,time(ct%),uprime(0,ct%),vprime(0,ct%)
  614.         'lprint "70",ct%,"time(ct%)="time(ct%),"uprime(0,ct%)="uprime(0,ct%)
  615.         'lprint "vprime(0,ct%)=",vprime(0,ct%):delay .3
  616.         ct% = ct% + 1                  'next time based data point
  617.         if ct% > 500 then goto 66
  618.         tst = tst + adt                 'next time step
  619.         next i%
  620.         goto 65
  621. 71      'write #2,ct%,time(ct%),uprime(0,ct%),vprime(0,ct%)
  622.         ct% = ct% + 1
  623.         if ct% > 500 then goto 66
  624.         tst = tst + adt
  625. 65 next k%
  626. 66
  627.  
  628.             ' uprime(x,ct%) => position at x (0=>surface) for ct% = 1 to 500
  629.             ' vprime(x,ct%) => strain at x (0 =>surface) for ct% = 1 to 500
  630.  
  631.      ' ... Obtain polished rod velocity by differentiating position data ...
  632.  
  633.    print " ... Data loaded and converted to time and strain ...
  634.  
  635.    ct% = ct% - 1
  636.  
  637.    for i% = 1 to ct%
  638.         if i%  = ct% then uprime(0,i%+1) = uprime(0,1)
  639.         ww(i%) = (uprime(0,i%+1) - uprime(0,i%))/adt   ' delta x / delta t
  640.    next i%
  641.  
  642.    wprime(0,1) = (ww(1) + ww(ct%))/2!
  643.  
  644.    for ii% = 2 to ct%
  645.         wprime(0,ii%) = (ww(ii%) + ww(ii%-1))/2!
  646.         'print ii%,time(ii%),wprime(0,ii%):delay .3
  647.    next ii%
  648.    'print "ct = 499"
  649.    'print uprime(0,499),vprime(0,499),wprime(0,499)
  650.    'print "ct = 500"
  651.    'print uprime(0,500),vprime(0,500),wprime(0,500)
  652.    'print "ct = 1"
  653.    'print uprime(0,1)vprime(0,1)wprime(0,1)
  654.  
  655.    '   ... Start method of characteristics ...
  656.    cls
  657.    print "Executing method of characteristics"
  658.  
  659. 75 for ijk% = 1 to ntap%
  660.  
  661.         ' Calculate the propagation interval dx(inches), number of propagation
  662.         ' intervals nn%, and fractional length for interpolation aaleng%(in),
  663.         ' for current rod section ijk%
  664.  
  665.         'print "adt=",adt,"a(ijk%)=",a(ijk%),"ijk%",ijk%
  666.  
  667.         dx = adt * a(ijk%) / 24! ':print dx
  668.         dx = dx * 12!
  669.         'print aleng(ijk%)
  670.         anum = aleng(ijk%) / dx       'fractional # of dx's in this section
  671.         nn% = int(anum)               'number of dx's in this section
  672.         ill% = 1
  673.         aaleng(ijk%) = aleng(ijk%) - (dx * nn%)   'remaining section past
  674.                                                   'nn%
  675.         'print : print anum,nn%,aaleng(ijk%)
  676.  
  677.     aa = ( 1! / a(ijk%)) + (( c * adt/2!)/(2! * a(ijk%)))
  678.     bb = ( 1! / a(ijk%)) - (( c * adt/2!)/(2! * a(ijk%)))
  679.  
  680.     '   ... Set convergence counter to zero ...
  681.  
  682.     lll% = 0
  683.  
  684.     '   ... Check if number of propagation intervals is zero ...
  685.     '       If so, go to interpolation
  686.  
  687.     if nn% = 0 then goto 76 else goto 80
  688.  
  689.  
  690. 76  for i% = 1 to ct%
  691.          uprime(ijk%,i%) = uprime(ijk%-1,i%)       'initial position
  692.          vprime(ijk%,i%) = vprime(ijk%-1,i%)       'initial strain
  693.          wprime(ijk%,i%) = wprime(ijk%-1,i%)       'initial velocity
  694.          'print "76",uprime(ijk%,i%),vprime(ijk%,i%),wprime(ijk%,i%)
  695.     next i%
  696. '   delay 20
  697.     goto 200
  698.  
  699. 80  for j% = 1 to ct%
  700.          if j% < ct% goto 90
  701.               uprime(ijk%-1,ct%+1) = uprime(ijk%-1,1)       'position
  702.               vprime(ijk%-1,ct%+1) = vprime(ijk%-1,1)       'strain
  703.               wprime(ijk%-1,ct%+1) = wprime(ijk%-1,1)       'velocity
  704.  
  705. 90       cc = - bb * wprime(ijk%-1,j%) + vprime(ijk%-1,j%)
  706.          dd =   aa * wprime(ijk%-1,j%+1) + vprime(ijk%-1,j%+1)
  707.  
  708.       '   ... calculate velocity and strain at next point ...
  709.  
  710.            wprime(ijk%,j%) = (dd - cc) / (bb + aa)
  711.            vprime(ijk%,j%) = (aa*dd + bb*cc) / (bb + aa)
  712.  
  713.       '   ...calculate positions at next two points along the char. ...
  714.  
  715. 92         u=(dx/2!)*(vprime(ijk%,j%)+vprime(ijk%-1,j%+1))+uprime(ijk%-1,j%+1)
  716.            u=u-.5*(adt/2!)*(wprime(ijk%,j%)+wprime(ijk%-1,j%+1))
  717.            uu=(dx/2!)*(vprime(ijk%,j%)+vprime(ijk%-1,j%))+uprime(ijk%-1,j%)
  718.            uu=uu+.5*(adt/2!)*(wprime(ijk%,j%)+wprime(ijk%-1,j%))
  719.            diff = u - uu
  720.            diff = abs(diff)
  721.  
  722.      'Check convergence of position along both characteristics. If convergence
  723.      'is satisfied, go on.  If not increase counter by 1 and iterate
  724.  
  725. 94         if diff < res goto 100
  726.            lll% = lll% + 1         'increase counter
  727.            f = ((c*adt)/(4!*a(ijk%)))*(wprime(ijk%,j%)+wprime(ijk%-1,j%))
  728.            f = f - wprime(ijk%-1,j%)/a(ijk%)+vprime(ijk%-1,j%)
  729.            g = ((c*adt)/(4!*a(ijk%)))*(wprime(ijk%,j%)+wprime(ijk%-1,j%+1))
  730.            g = g + wprime(ijk%-1,j%+1)/a(ijk%)+vprime(ijk%-1,j%+1)
  731.  
  732.            vprime(ijk%,j%) = (f+g)/2!
  733.            wprime(ijk%,j%) = a(ijk%) * (g-f)/2!
  734.  
  735.                 'if more than ten iterations are needed to converge, prompt
  736.                 'for larger value of res(inches).
  737.            if lll% > 10 then goto 95 else goto 92
  738.  
  739. 95              print "Convergence criteria, RES, too small."
  740.                 print "Currently RES in inches is",res
  741.                 input "Enter new value for RES:",res
  742.                      goto 75
  743.  
  744. 100     lll% = 0
  745.         uprime(ijk%,j%) = (u+uu)/2
  746.         'if j% >498 then print uprime(ijk%,j%),vprime(ijk%,j%),wprime(ijk%,j%)
  747.     next j%
  748.     'end of loop to calculate all points (ct%) at current dx station
  749.  
  750.     if ill% < nn% then goto 110 else goto 200
  751.  
  752. 110      ill% = ill%+1
  753.          for i% = 1 to ct%
  754.               uprime(ijk%-1,i%) = uprime(ijk%,i%)    'position
  755.               vprime(ijk%-1,i%) = vprime(ijk%,i%)    'strain
  756.               wprime(ijk%-1,i%) = wprime(ijk%,i%)    'velocity
  757.          next i%
  758.      goto 80
  759.  
  760.      'End of loop to calculate all data at nn% propagation intervals for
  761.      'current rod section.
  762.  
  763. 200  'check if interpolation is desired on current rod section ijk%
  764.      if aaleng(ijk%) < cut goto 300
  765.      ddt = aaleng(ijk%) / a(ijk%)
  766.      dt = adt
  767.  
  768.      us(1) =(1!/dt)*(uprime(ijk%,ct%)*ddt+uprime(ijk%,1)*(dt-ddt))
  769.      vs(1) =(1!/dt)*(vprime(ijk%,ct%)*ddt+vprime(ijk%,1)*(dt-ddt))
  770.      ws(1) =(1!/dt)*(wprime(ijk%,ct%)*ddt+wprime(ijk%,1)*(dt-ddt))
  771.  
  772.      lm% = 2 * ct%
  773.  
  774.      us(lm%) = (1!/dt)*(uprime(ijk%,ct%)*(dt-ddt)+uprime(ijk%,1)*ddt)
  775.      vs(lm%) = (1!/dt)*(vprime(ijk%,ct%)*(dt-ddt)+vprime(ijk%,1)*ddt)
  776.      ws(lm%) = (1!/dt)*(wprime(ijk%,ct%)*(dt-ddt)+wprime(ijk%,1)*ddt)
  777.  
  778.      lc% = lm% - 2
  779.      ic% = 0
  780.      for k% = 2 to lc% step 2
  781.           iaa% = k%-1-ic%
  782.           ia%  = k%-ic%
  783.  
  784.           us(k%) = (1!/dt)*(uprime(ijk%,iaa%)*(dt-ddt)+uprime(ijk%,ia%)*ddt)
  785.           vs(k%) = (1!/dt)*(vprime(ijk%,iaa%)*(dt-ddt)+vprime(ijk%,ia%)*ddt)
  786.           ws(k%) = (1!/dt)*(wprime(ijk%,iaa%)*(dt-ddt)+wprime(ijk%,ia%)*ddt)
  787.  
  788.           us(k%+1)=(1!/dt)*(uprime(ijk%,iaa%)*ddt+uprime(ijk%,ia%)*(dt-ddt))
  789.           vs(k%+1)=(1!/dt)*(vprime(ijk%,iaa%)*ddt+vprime(ijk%,ia%)*(dt-ddt))
  790.           ws(k%+1)=(1!/dt)*(wprime(ijk%,iaa%)*ddt+wprime(ijk%,ia%)*(dt-ddt))
  791.  
  792.           ic% = ic% + 1
  793.     next k%
  794.  
  795.     'for i% = 1 to lm%
  796.     '     print us(i%),vs(i%),ws(i%)
  797.     'next i%
  798.  
  799.     aaa = (1!/a(ijk%)) + ((c*ddt)/( 2! * a(ijk%)))
  800.     bbb = (1!/a(ijk%)) - ((c*ddt)/( 2! * a(ijk%)))
  801.  
  802.     for k% = 1 to ct%
  803.  
  804.          kk% = 2 * k%
  805.          kkk% = kk% - 1
  806.  
  807.          ccc = -bbb * ws(kkk%) + vs(kkk%)
  808.          ddd = +aaa * ws(kk%) + vs(kk%)
  809.  
  810.          'If this is the last section, the following calculated positions
  811.          'velocities, and strains will be the conditions at the pump.
  812.  
  813.          wprime(ijk%,k%) = (ddd - ccc)/(bbb + aaa)
  814.          vprime(ijk%,k%) = (aaa * ddd + bbb * ccc)/(bbb + aaa)
  815.  
  816.          uans = (aaleng(ijk%)/2!)*(vprime(ijk%,k%)+vs(kk%))
  817.          uans = uans + us(kk%) - (ddt*(wprime(ijk%,k%) + ws(kk%))/2!)
  818.          uuans = (aaleng(ijk%)/2!) * (vprime(ijk%,k%) + vs(kkk%))
  819.          uuans = uuans + us(kkk%) + (ddt*(wprime(ijk%,k%) + ws(kkk%))/2!)
  820.          uprime(ijk%,k%) = (uans + uuans)/2!
  821.  
  822.          next k%
  823.  
  824.     '   Check to see if this is the last rod section. If not use continuity
  825.     '   to determine boundary conditions.   # of sections loop
  826.     '   ends at 310
  827.  
  828. 300  if ijk% < ntap% then goto 305 else goto 310
  829. 305  for i% = 1 to ct%
  830.  
  831.           vprime(ijk%,i%) = vprime(ijk%,i%) * slfact(ijk%)  'load in this sect
  832.           vprime(ijk%,i%) = vprime(ijk%,i%) / slfact(ijk%+1) 'continuity load
  833.           'print i%
  834.           'print uprime(ijk%,i%),vprime(ijk%,i%),wprime(ijk%,i%) 'converted to
  835.                                                              'strain in next
  836.                                                              'section
  837.      next i%
  838.  
  839. 310  next ijk%      'finish method of characteristics
  840.  
  841.      for i% = 1 to ct%
  842.  
  843. 330        uprime(ntap%,i%) = -uprime(ntap%,i%) + rodst
  844.            vprime(ntap%,i%) = vprime(ntap%,i%)  * slfact(ntap%)
  845.  
  846.      next i%
  847.  
  848.       pstmin   =   1000
  849.       pstmax   =  -1000
  850.       loadmax  =   0
  851.       effcount =  0
  852.       loadmin  = 10000000
  853.  
  854.       for i% = 1 to ct%
  855.            if uprime(ntap%,i%) < pstmin then pstmin = uprime(ntap%,i%)
  856.            if uprime(ntap%,i%) > pstmax then pstmax = uprime(ntap%,i%)
  857.            if vprime(ntap%,i%) > loadmax then loadmax = vprime(ntap%,i%)
  858.            if vprime(ntap%,i%) < loadmin then loadmin = vprime(ntap%,i%)
  859.  
  860.            if i% > (ct%/2) then goto 341 else goto 343
  861. 341          if effcount <> 0 then goto 343
  862.              if vprime(ntap%,i%) < .05 * loadmax then goto 342 else goto 343
  863. 342          effstr = uprime(ntap%,i%)
  864.              effcount = 1
  865. 343   next i%
  866.  
  867.       lprint"     Downhole Dynagraph Data"
  868.       lprint using"          Maximum Pumpstroke (inches)=###";pstmax - pstmin
  869.       lprint
  870.  
  871.       effstr = effstr - pstmin
  872.       lprint using"          Effective Pumpstroke (inches) =###";effstr:lprint
  873.  
  874.       thpumprt =  effstr * nst# * plunarea * 60 * 24/(1728*5.615)
  875.       lprint using"          Theoretical Pumping Rate (BBL/day)=####";thpumprt
  876.       lprint
  877.  
  878.       s = 0
  879.       peaktk = 0:mintk =1000000
  880.       dim torque(500)
  881.       dim angle(500)
  882.       pktkang = -1
  883.       for i% = 1 to 500
  884.            if i% > 250 then s = 1
  885.  
  886.            sine   = posi(i%)/halfstr
  887.            cosine = ((1-(sine)^2)^.5)
  888.  
  889.            if (cosine = 0 and posi(j%) < 0) then angle(i%) = -90  'avoid
  890.            if (cosine = 0 and posi(j%) > 0) then angle(i%) =  90  'error
  891.  
  892.            if cosine <> 0 then angle(i%) = (atn((sine)/(cosine)))*180 /pi
  893.            if s = 0 then angle(i%) = angle(i%) + 90
  894.            if s = 1 then angle(i%) = 270 - angle(i%)
  895.            If i% = 1 then angle(i%) = 0
  896.            If i% = 500 then angle(i%) = 360
  897.            torque(i%) = ( load(i%) - cbe% ) * halfstr * sin(angle(i%)*2*pi/360)
  898.            if torque(i%) > peaktk then pktkang = angle(i%)
  899.            if torque(i%) > peaktk then peaktk  = torque(i%)
  900.            if torque(i%) < mintk then mintk  = torque(i%)
  901.            'print torque(i%),angle(i%):delay .1
  902.       next i%
  903.  
  904.      lprint "     *********************************":lprint
  905.      lprint "     Torque Data"
  906.      lprint using"          Peak Torque (inch lb) = ########"; peaktk :lprint
  907.      lprint using"          Peak Torque Angle (degrees) =###"; pktkang :lprint
  908.      lprint "     *********************************":lprint
  909.  
  910.  
  911. input "Do you have a pen plotter available";ans$
  912.       if ans$ = "y" then goto plot:
  913.       if ans$ = "Y" then goto plot:
  914. cls
  915.  
  916. print "Your downhole dynagraph data will be sent to a LOTUS 1-2-3 file"
  917. print "where a graph can be made."
  918.  
  919. 340   chdir"c:\123\files"
  920.       files "*.prn"
  921.  
  922.       input"LOTUS file for Downhole Pump Card?",file$
  923.  
  924.       open file$ for output as #2
  925.            for i% = 1 to ct%
  926.                 write #2,i%,time(i%),uprime(ntap%,i%),vprime(ntap%,i%),angle(i%),torque(i%)
  927.            next i%
  928.       close #2
  929.       lprint "      Data saved under:",file$
  930. input "Press enter to return to system.",ent
  931.       cls
  932.       end
  933.  
  934. plot:
  935.      print"You may need to switch from the digitizer to the plotter."
  936.      print"The plotter will need pens in both the 1st and 2cd slots."
  937.      input"Press enter when ready."ans$
  938.  
  939.      loadmax = ceil(loadmax/500) * 500   'set max load as nearest 500 above
  940.      loadmin = int(loadmin/500) * 500    'set min load as nearest 500 below
  941.  
  942.      pstmax  = ceil(pstmax/5) * 5        'set max position as nearest 5 above
  943.      pstmin  = int(pstmin/5) * 5         'set min position as nearest 5 below
  944. close #1
  945.  
  946. open "com1:9600 ,n,8,1,rs,cs65535,ds,cd" as #1
  947. print #1, "in;sp1;ip2000,2000,8750,6250;"
  948. print #1, "sc";pstmin,pstmax,loadmin,loadmax;
  949. print #1, "pu"pstmin,loadmin;"pd"pstmax,loadmin,pstmax,loadmax,pstmin,loadmax,pstmin,loadmin;"pu;"
  950. print #1, "si.2,.3;tl1.5,0;"
  951.  
  952. for x = pstmin to pstmax step 10
  953.       print #1,"pa";x,loadmin;";xt;"
  954.       if x > 100 then print #1,"cp -2.27,-1;":goto 800
  955.       if x < 10 then print #1,"cp -1.33,-1; "
  956.       if x > 9 then print #1,"cp -1.80,-1;"
  957. 800      print #1,"lb";x;chr$(3)
  958.  
  959. next x
  960.  
  961. pstavg = (pstmax - pstmin)/2 + pstmin
  962. print #1,"pa"pstavg,loadmin;"cp-9,-2.5;"
  963. print #1,"lbPump Position - in"+chr$(3)
  964.  
  965. for x = loadmin to loadmax step 500
  966.       print #1,"pa";pstmin,x;";yt;"
  967. next x
  968.  
  969. for x = loadmin to loadmax step 1000
  970.   print #1,"pa"pstmin,x;
  971.      if x < -9999 then print #1,"cp-7,-.25;"
  972.      if x <  -999 then print #1,"cp-6,-.25;":goto 801
  973.      if x <  -100 then print #1,"cp-5,-.25;":goto 801
  974.      if x =     0 then print #1,"cp-2,-.25;":goto 801
  975.      if x =   500 then print #1,"cp-4,-.25;":goto 801
  976.      if x >  9999 then print #1,"cp-6,-.25;":goto 801
  977.      if x >   999 then print #1,"cp-5,-.25;":goto 801
  978.      if x > 99999 then print #1,"cp-7,-.25;":goto 801
  979.  
  980. 801  print #1,"lb"x;chr$(3)
  981. next x
  982.  
  983. print #1,"pa"pstmin,loadmax;"cp-5,2;"
  984. print #1,"lbLoad - lb"+Chr$(3)"
  985.  
  986. title$ = wellid$
  987.  
  988. print #1,"sp1;pa"pstavg,loadmax;"si.3,.5;cp"; -len(title$)/2;"2;lb";title$+chr$(3)
  989.  
  990.  
  991. print #1,"sp1;si.2,.3;"
  992.  
  993. print #1,"sp2;lt;pa";uprime(ntap%,1)","vprime(ntap%,1);"pd";
  994.      for x = 1 to 500
  995.            print #1,"pd";uprime(ntap%,x)","vprime(ntap%,x);
  996.      next x
  997.  
  998.  
  999. 6101  print #1,"sp0;"
  1000. cls
  1001.  
  1002. print"Reload plotter paper to plot Torque Curve."
  1003. input"Press enter when ready."ans$
  1004.  
  1005.      loadmax = ceil(peaktk/5000) * 5   'set max torque as nearest 5000 above
  1006.      loadmin = int(mintk/5000) * 5    'set min torque as nearest 5000 below
  1007.      if loadmin > 0 then loadmin = 0
  1008.  
  1009.      pstmax  = 360
  1010.      pstmin  = 0
  1011. close #1
  1012. open "com1:9600 ,n,8,1,rs,cs65535,ds,cd" as #1
  1013. print #1, "in;sp1;ip2000,2000,8750,6250;"
  1014. print #1, "sc";pstmin,pstmax,loadmin,loadmax;
  1015. print #1, "pu"pstmin,loadmin;"pd"pstmax,loadmin,pstmax,loadmax,pstmin,loadmax,pstmin,loadmin;"pu;"
  1016. print #1, "si.2,.3;tl1.5,0;"
  1017. for x = pstmin to pstmax step 45
  1018.       print #1,"pa";x,loadmin;";xt;"
  1019.       if x > 100 then print #1,"cp -2.27,-1;":goto 805
  1020.       if x < 10  then print #1,"cp -1.33,-1; "
  1021.       if x > 9   then print #1,"cp -1.80,-1;"
  1022. 805      print #1,"lb";x;chr$(3)
  1023. next x
  1024. pstavg = 180
  1025. print #1,"pa"pstavg,loadmin;"cp-7,-2.5;"
  1026. print #1,"lbAngle - degrees"+chr$(3)
  1027.  
  1028. for x = loadmin to loadmax step 2.5
  1029.       print #1,"pa";pstmin,x;";yt;"
  1030. next x
  1031. for x = loadmin to loadmax step 5
  1032.   print #1,"pa"pstmin,x;
  1033.      if x < -9.999 then print #1,"cp-3,-.25;":goto 806
  1034.      if x <  -.999 then print #1,"cp-2,-.25;":goto 806
  1035.      if x <  -.100 then print #1,"cp-2,-.25;":goto 806
  1036.      if x =     0 then print #1,"cp-2,-.25;":goto 806
  1037.      if x =   .500 then print #1,"cp-2,-.25;":goto 806
  1038.      if x > 99.999 then print #1,"cp-4,-.25;":goto 806
  1039.      if x >  9.999 then print #1,"cp-3,-.25;":goto 806
  1040.      if x >   .999 then print #1,"cp-2,-.25;":goto 806
  1041.  
  1042.  
  1043. 806  print #1,"lb"x;chr$(3)
  1044. next x
  1045.  
  1046. print #1,"pa"pstmin,loadmax;"cp-7.5,2;"
  1047. print #1,"lbTorque - ( M in-lb)"+Chr$(3)"
  1048.  
  1049. title$ = wellid$
  1050.  
  1051. print #1,"sp1;pa"pstavg,loadmax;"si.3,.5;cp"; -len(title$)/2;"2;lb";title$+chr$(3)
  1052.  
  1053.  
  1054. print #1,"sp1;si.2,.3;"
  1055.  
  1056. print #1,"sp2;lt;pa";angle(1)","torque(1);"pd";
  1057.      for x = 1 to 500
  1058.            print #1,"pd";angle(x)","torque(x)/1000;
  1059.      next x
  1060.  
  1061.  
  1062.  print #1,"sp0;"
  1063.  
  1064. chdir "c:\turbo\files"
  1065. close #1
  1066. end
  1067.