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