FullW # 1
FontName = "Courier"
FontSize = 10
PrintScroll = 1
AutoRedraw = 1
Global Int nz, np, i, nss, ncal, dd, yd, n, n0, np1, np2, xx, yy
Global Int nz1, nz2, iflog, r, rn, icase, nlfp, nz0, zc
Global Double a0, fnorm, zre(99), zim(99), pre(99), pim(99)
Global Double sens, fsens, dmp, dn, g, om, per, sq, sw, trans
Global Double zq, dq, f, flog, freq, ymax, pi2, pih, pifac, y
Global Double s(99), sre(99), sim(99), s0(99), s0re(99), s0im(99)
Global Double ph, yph, yph1, yph2, xxa, yya, ypha, yph1a, yph2a
Global Double yph3, yph3a, slope, slopa, floga, pha, flattest
Global Double tnorm, a0poles, tsens, a0sens, sqmin, lfpre, lfpim
Global Double lfpq, tas, a1, atn0, phi0, omn, t1hz
Global String fname, zeile, code, par, ok, sta, hz, typ, outfile
Global String ig, ex, sc, unitsin, unitsout, net, pp, cdir
Global Boolean vresp, dresp

cdir = CurDir$
Print
Print " WINRESP evaluates a SEED-format RESP file (the text version of headerless SEED)"
Print " for the response of the sensor ('stage 1'). It produces a graph of the amplitude"
Print " and phase response on the screen, a file with the numerical values of the response"
Print " at frequencies of your choice, and will help to determine the correct entries for"
Print " the normalization factor and the sensitivity."
Print " Usage:"
Print " Start from an existing RESP file. It must be in a subdirectory RESPfiles of your"
Print " working directory. Edit zeros, poles, and sensitivity in the RESP file as far as"
Print " you know their values from a calibration or from a manual. If you don't know what"
Print " to enter as the sensitivity of a broadband seismometer, enter the nominal generator"
Print " constant although this is in general not the correct entry here. The 'frequency of"
Print " calibration' should be in the middle of the passband and agree with the normalization"
Print " frequency, but this is not essential if the sensitivity parameter is correct."
Print
Print " Prepare a list of frequencies, one per line, with the name 'frequencies.lst' in the"
Print " working directory if you need numerical values of the response. They will be stored"
Print " in a file 'response.lst' in the working directory."
Print
Print " A screen-dump program like Gadwin PrintScreen will be useful to save the response graph."
Print
Input "=============> (1) ok?   ";ok
Repeat
  Print
  FileSelect "RESPfiles\*", "", fname
  If fname = ""
    Print "no file selected"
    Stop
  EndIf
  ' read RESP file
  zc = 0
  Open fname for Input As # 2
  Repeat
    Inc zc
    If zc > 30 Then
      zc = 0
      Print
      Input "=============> (2) ok?   ";ok
    EndIf
    Line Input # 2, zeile
    Print zeile
    code = Left(zeile, 7)
    If code = "B050F03"
      sta = Mid(zeile, 25)
    Else If code = "B050F16"
      net = Mid(zeile, 25)
    Else If code = "B053F03" Or code = "B054F03"
      typ = Mid(zeile, 51, 3)
    Else If code = "B053F05"
      unitsin = Mid(zeile, 51)
    Else If code = "B053F06"
      unitsout = Mid(zeile, 51)
    Else If code = "B053F07"
      a0 = Val(Mid(zeile, 51))
    Else If code = "B053F08"
      fnorm = Val(Mid(zeile, 51))
    Else If code = "B053F09"
      nz = Val(Mid(zeile, 51))
    Else If code = "B053F14"
      np = Val(Mid(zeile, 51))
    Else If code = "B053F10"
      par = Mid(zeile, 12)
      i = Val(par)
      par = Mid(par, Val?(par) + 1)
      zre(i + 1) = Val(par)
      par = Mid(par, Val?(par) + 1)
      zim(i + 1) = Val(par)
    Else If code = "B053F15"
      par = Mid(zeile, 12)
      i = Val(par)
      par = Mid(par, Val?(par) + 1)
      pre(i + 1) = Val(par)
      par = Mid(par, Val?(par) + 1)
      pim(i + 1) = Val(par)
    Else If code = "B058F03"
      nss = Val(Mid(zeile, 51))
    Else If code = "B058F04"
      sens = Val(Mid(zeile, 51))
    Else If code = "B058F05"
      fsens = Val(Mid(zeile, 51))
    Else If code = "B058F06"
      ncal = Val(Mid(zeile, 51))
    EndIf
  Until code = "B058F06" Or EOF(# 2)

  Print
  a0 = Abs(a0)
  fnorm = Abs(fnorm)
  fsens = Abs(fsens)
  Input "=============> (3) ok?   ";ok
  Print " parameters read:"
  Print " a0,fnorm,nz,np: "; a0; fnorm; nz; np
  If InStr(typ, "A") > 0
    hz = "r"
  Else If InStr(typ, "B") > 0
    hz = "h"
  Else
    Print " ==> unknown response type "; typ; ", can handle only A and B"
    Stop
  EndIf
  pi2 = 2. * PI
  pih = PI / 2.
  If hz = "r"
    pifac = pi2
    Print " poles and zeros given as angular frequencies (rad/s)"
  Else
    pifac = 1.
    Print " poles and zeros given as common frequencies (Hz)"
  EndIf
  Print " Zeros (re and im parts):"
  For i = 1 To nz
    Print " "; zre(i), zim(i)
  Next i
  Print " Poles (re and im parts):"
  For i = 1 To np
    Print " "; pre(i), pim(i)
    If pre(i) >= 0
      Print
      Print " Error in the RESP file. Poles must lie in the left halfplane (Re s<0). STOP."
      Stop
    EndIf
  Next i
  Print " nss,sens,fsens,ncal:   "; nss; sens; fsens; ncal
  If sens > 0
    phi0 = 0
  Else
    sens = -sens
    phi0 = PI
    Print " Negative sensitivity will be expressed by the phase"
  EndIf
  Close # 2
  Print
  ' interface to POLZERO input format
  ' =================================
  ' identify single (real) and paired (complex-conjugate) zeros
  ' single zeros are stored as s0(i), i=1..nz1
  ' paired zeros are stored as s0re(i), s0im(i), i=1..nz2
  nz1 = 0
  nz2 = 0
  For i = 1 To nz
    If zre(i + 1) = zre(i) And zim(i) <> 0 And zim(i + 1) = -zim(i)
      Inc nz2
      Inc i
      s0re(nz2) = zre(i)
      s0im(nz2) = zim(i)
    Else
      If zim(i) <> 0
        Print "A single zero must be real: stop"
        Stop
      EndIf
      Inc nz1
      s0(nz1) = zre(i)
    EndIf
  Next i
  ' identify single (real) and paired (complex-conjugate) poles
  ' single poles are stored as s(i), i=1..np1
  ' paired poles are stored as sre(i), s0im(i), i=1..np2
  np1 = 0
  np2 = 0
  sqmin = 1e12
  For i = 1 To np
    If pre(i + 1) = pre(i) And pim(i + 1) = -pim(i)
      Inc np2
      Inc i
      sre(np2) = pre(i)
      sim(np2) = pim(i)
      If pre(i) ^ 2 + pim(i) ^ 2 < sqmin
        sqmin = pre(i) ^ 2 + pim(i) ^ 2
        nlfp = np2
        lfpre = pre(i)
        lfpim = pim(i)
      EndIf
    Else
      If pim(i) <> 0
        Print "A single pole must be real: STOP"
        Stop
      EndIf
      Inc np1
      s(np1) = pre(i)
    EndIf
  Next i
  Input "=============> (4) ok?   ";ok

  ' generalized POLZERO routine
  ' ===========================
  Print
  Print " Interpreting RESP-file: "; fname
  Print " ===================================="
  nz0 = 0
  If nz1 > 0
    Print
    Print " Single Zeros:"
    For n = 1 To nz1
      Print " Zero #"; n; ": "
      freq = -s0(n) / pifac
      Print " s", Spc(1), "freq", "per"
      If freq = 0
        Inc nz0
        Print s0(n), Spc(1), Round(freq, 3)
      Else
        per = 1 / freq
        Print s0(n), Spc(1), Round(freq, 3), Round(per, 3)
        If freq < 0
          Print " >>>>> system is mixed-delay - is this OK or a typing error?"
        EndIf
      EndIf
    Next n
  EndIf
  If nz2 > 0
    Print " Paired Zeros:"
    For n = 1 To nz2
      Print " Pair #"; n; ": "
      sq = s0re(n) ^ 2 + s0im(n) ^ 2
      sw = Sqr(sq)
      freq = sw / pifac
      per = 1 / freq
      dmp = -sre(n) / sw
      Print " sre", "+-sim", "freq", "per", "dmp"
      Print s0re(n), -s0im(n), Round(freq, 3), Round(per, 3), Round(dmp, 4)
    Next n
  EndIf
  If np1 > 0
    Print
    Input "============= (5) ok?   ";ok
    Print " Single Poles:"
    For n = 1 To np1
      Print " Pole #"; n; ": "
      Exit If s(n) = 0
      freq = -s(n) / pifac
      per = 1 / freq
      Print " s1", Spc(1), "freq", "per"
      Print s(n), Spc(1), Round(freq, 3), Round(per, 3)
    Next n
  EndIf
  If np2 > 0
    Print
    Input "============= (6) ok?   ";ok
    Print " Pole Pairs:"
    For n = 1 To np2
      Print " Pair #"; n; ": "
      sq = sre(n) ^ 2 + sim(n) ^ 2
      sw = Sqr(sq)
      freq = sw / pifac
      per = 1 / freq
      dmp = -sre(n) / sw
      Print " sre", "+-sim", "freq", "per", "dmp"
      Print " "; sre(n), -sim(n), Round(freq, 3), Round(per, 3), Round(dmp, 4)
    Next n
    Print
  EndIf
  Print
  vresp = (nz0 = 2) And np2 > 0
  If vresp
    Print " Pole pair #"; nlfp; "together with the two zeros at zero frequency is assumed"
    Print " to represent the high-pass section of a flat response to velocity."
    lfpq = lfpre ^ 2 + lfpim ^ 2
  EndIf
  dresp = (nz0 = 3) And np2 > 0
  Print
  Input "=============> (7) ok?   ";ok
  ' calculate the response at the normalization frequency, at the freq. of sensitivity,
  ' and the asymptotic response
  For icase = 1 To 3
    If icase = 1
      om = Max(pifac * fnorm, 1e-24)
      omn = om
    Else If icase = 2
      om = Max(pifac * fsens, 1e-24)
    Else
      om = pifac
    EndIf
    zq = 1
    For n = 1 To nz1
      zq = zq * (s0(n) ^ 2 + om ^ 2)
    Next n
    For n = 1 To nz2
      zq = zq * (s0re(n) ^ 2 + (om + s0im(n)) ^ 2)
      zq = zq * (s0re(n) ^ 2 + (om - s0im(n)) ^ 2)
    Next n
    dq = 1
    For n = 1 To np1
      dq = dq * (s(n) ^ 2 + om ^ 2)
    Next n
    For n = 1 To np2
      dq = dq * (sre(n) ^ 2 + (om + sim(n)) ^ 2)
      dq = dq * (sre(n) ^ 2 + (om - sim(n)) ^ 2)
    Next n
    If icase = 1
      tnorm = Sqr(zq / dq)
      a0poles = 1. / tnorm
      tas = Sqr(((lfpq - om ^ 2) ^ 2 + 4. * om ^ 2 * lfpre ^ 2) / om ^ 4)
    Else If icase = 2
      tsens = Sqr(zq / dq)
    Else
      t1hz = Sqr(zq / dq)
    EndIf
  Next icase
  Print
  Print " Using the normalization factor A0 = "; pp(a0);
  Print  " from line B053F07 of the RESP file,"
  Print " the normalized response at"; fnorm; "Hz is "; pp(a0 * tnorm)
  Print
  Print " If this is not very close to 1 then A0 in the RESP file is wrong."
  a0 = a0poles
  Print " The exact value can be calculated from the poles and zeros and is ";
  Print " A0 = "; pp(a0)
  Print
  Print " The exact value will be used in the following, and should be entered into the"
  Print " RESP file to avoid ambiguity; evaluation programs may use one value or the other."
  Print
  tsens = a0 * tsens ' tnorm now equals 1 and no longer appears in the equations
  t1hz = a0 * t1hz
  Print " If the sensitivity in line B058F04,"; pp(sens); " is in fact the sensitivity"
  Print " measured at the 'frequency of sensitivity',"; fsens; "Hz, then: "
  Print " the gain at the normalization frequency,"; fnorm; "Hz, is"; pp(sens / tsens)
  If vresp
    Print " the asymptotic high-frequency gain of the high-pass section is";
    Print pp(sens / tsens * tas)
  EndIf
  Print
  Print " If the sensitivity in line B058F04 WAS NOT measured at the frequency of sensi-"
  Print " tivity, then the entry should be changed if the difference is substantial."
  Print " What kind of sensitivity is really given the existing RESP file?"
  Print " the gain at the frequency of Sensitivity "; fsens; "Hz  -     answer s"
  Print " the gain at the Normalization frequency "; fnorm; "Hz -       answer n"
  If vresp
    Print " the Asymptotic high-frequency gain of the high-pass section"
    Print " (the nominal generator constant is usually of this kind) -    answer a";
  EndIf
  Repeat
    Input " ? "; ig
    ig = Lower(ig)
  Until ig = "s" Or ig = "" Or ig = "n" Or (vresp And ig = "a")
  If ig = "s" Or ig = ""
    g = sens / tsens
  Else If ig = "n"
    g = sens
  Else
    g = sens / tas
  EndIf
  ' g is now the gain factor that must be multiplied with the normalized response.
  Print
  Print " The 'sensitivity' in line B058F04 of the RESP file must be set to"; pp(g * tsens)
  Print
  Print " This value will be used for the graph and for the file response.lst"
  Print
  Print " With the new values,"
  Print " the gain at the frequency of sensitivity,"; fsens; "Hz, is"; pp(g * tsens)
  Print " the gain at the normalization frequency,"; fnorm; "Hz, is"; pp(g)
  If vresp
    Print " the asymptotic high-frequency gain of the high-pass section is";
    Print pp(g * tas)
  EndIf
  Print
  Print " =========================================================================="
  Print "                                SUMMARY"
  Print " The normalization factor should be "; pp(a0); " and the sensitivity"; pp(g * tsens)
  Print " =========================================================================="
  Print
  If dresp
    Print " In order to convert the response to ground displacement into the response to"
    Print " ground velocity, the RESP file should be modified as follows:"
    Print " the input units are m/s (line 'Response in units lookup')"
    Print " one zero at zero frequency must be cancelled"
    Print " the normalization factor A0 must be changed to "; pp(a0 * omn)
    Print " the sensitivity must be changed to "; pp(g * tsens / omn)
    Print " =========================================================================="
    Print
  EndIf

  Print " If a list of frequencies (one per line) is supplied as file frequencies.lst,"
  Print " then the response at these frequencies will be stored in file response.lst"
  Print
  Input "=============> (8) ok?   ";ok

  ' plot the response from 0.0001 Hz to 10000 Hz
  Cls
  ymax = -77777
  For iflog = 0 To 400
    flog = -4. + iflog / 50.
    f = 10. ^ flog
    om = pifac * f
    zq = 1.
    For n = 1 To nz1
      zq = zq * (s0(n) ^ 2 + om ^ 2)
    Next n
    For n = 1 To nz2
      zq = zq * (s0re(n) ^ 2 + (om + s0im(n)) ^ 2)
      zq = zq * (s0re(n) ^ 2 + (om - s0im(n)) ^ 2)
    Next n
    dq = 1.
    For n = 1 To np1
      dq = dq * (s(n) ^ 2 + om ^ 2)
    Next n
    For n = 1 To np2
      dq = dq * (sre(n) ^ 2 + (om + sim(n)) ^ 2)
      dq = dq * (sre(n) ^ 2 + (om - sim(n)) ^ 2)
    Next n
    ymax = Max(ymax, Log10(a0 * g * Sqr(zq / dq)))
  Next iflog
  dn = Round(5.45 - ymax)
  DefLine 0, 1
  QBColor 0
  Text 136, 660, "log f"
  QBColor 12
  Text 80, 385, "log"
  Text 80, 405, " R"
  QBColor 9
  Text 925, 385, "Phase"
  Text 925, 405, " *PI"
  QBColor 0
  For i = 20 To 820 Step 100
    Line i + 100, 55, i + 100, 655
    Text i + 83, 660, Str((i - 420) / 100)
  Next i
  For i = 10 To 610 Step 100
    Line 120, i + 45, 920, i + 45
    Text 90, i + 34, Str(Round((610 - i) / 100 - dn))
    Text 925, i + 34, Str(Round((610 - i) / 100) - 6)
  Next i
  ymax = -77777
  For iflog = 0 To 800
    flog = -4. + iflog / 100.
    f = 10. ^ flog
    om = pifac * f
    zq = 1.
    ph = phi0
    For n = 1 To nz1
      ph = ph - FN atn0(om, s0(n))
      zq = zq * (s0(n) ^ 2 + om ^ 2)
    Next n
    For n = 1 To nz2
      ph = ph - FN atn0(om + s0im(n), s0re(n))
      ph = ph - FN atn0(om - s0im(n), s0re(n))
      zq = zq * (s0re(n) ^ 2 + (om + s0im(n)) ^ 2)
      zq = zq * (s0re(n) ^ 2 + (om - s0im(n)) ^ 2)
    Next n
    dq = 1.
    For n = 1 To np1
      ph = ph + FN atn0(om, s(n))
      dq = dq * (s(n) ^ 2 + om ^ 2)
    Next n
    For n = 1 To np2
      ph = ph + FN atn0(om + sim(n), sre(n))
      ph = ph + FN atn0(om - sim(n), sre(n))
      dq = dq * (sre(n) ^ 2 + (om + sim(n)) ^ 2)
      dq = dq * (sre(n) ^ 2 + (om - sim(n)) ^ 2)
    Next n
    trans = a0 * g * Sqr(zq / dq)
    y = Log10(trans) + dn
    ymax = Max(ymax, y)
    ph = ph / PI
    xx = Round(520 + 100 * flog)
    yph = Round(55 - 100 * ph)
    yph1 = Min(655, Max(55, yph + 200))
    yph2 = Min(655, Max(55, yph + 400))
    yph3 = Min(655, Max(55, yph + 600))
    yph = Min(655, Max(55, yph))
    yy = Min(655, Max(55, Round(655 - 100 * y)))
    If iflog > 0 And (yy > 55 Or yya > 55) And (yy < 655 Or yya < 655)
      DefLine 0, 3
      QBColor 12
      Line xxa, yya, xx, yy
    EndIf
    If iflog > 0
      DefLine 0, 1
      QBColor 9
      Line xxa, ypha, xx, yph
      Line xxa, yph1a, xx, yph1
      Line xxa, yph2a, xx, yph2
      Line xxa, yph3a, xx, yph3
    EndIf
    xxa = xx
    yya = yy
    ypha = yph
    yph1a = yph1
    yph2a = yph2
    yph3a = yph3
    pha = ph
    floga = flog
    slopa = slope
  Next iflog
  QBColor 0
  Text 340, 25, "Network: " + net + "   Station: " + sta
  Text 120, 693, "IN: " + unitsin + ".  OUT: " + unitsout
  Text 120, 720, "gain at 1 Hz is " + pp(g * t1hz)
  Plot 460, 720
  Input " Type n for next or <cr> for end ", ok
Until ok = ""

' write response into file response.lst
' =====================================
ChDir cdir
If Exist("frequencies.lst")
  Open "frequencies.lst" for Input As # 3
  Open "response.lst" for Output As # 4
  Print # 4, "frequency", "amplitude", "phase in rad"
  While Not EOF(# 3)
    Input # 3, f
    om = pifac * f
    zq = 1.
    ph = phi0
    For n = 1 To nz1
      ph = ph - FN atn0(om, s0(n))
      zq = zq * (s0(n) ^ 2 + om ^ 2)
    Next n
    For n = 1 To nz2
      ph = ph - FN atn0(om + s0im(n), s0re(n))
      ph = ph - FN atn0(om - s0im(n), s0re(n))
      zq = zq * (s0re(n) ^ 2 + (om + s0im(n)) ^ 2)
      zq = zq * (s0re(n) ^ 2 + (om - s0im(n)) ^ 2)
    Next n
    dq = 1.
    For n = 1 To np1
      ph = ph + FN atn0(om, s(n))
      dq = dq * (s(n) ^ 2 + om ^ 2)
    Next n
    For n = 1 To np2
      ph = ph + FN atn0(om + sim(n), sre(n))
      ph = ph + FN atn0(om - sim(n), sre(n))
      dq = dq * (sre(n) ^ 2 + (om + sim(n)) ^ 2)
      dq = dq * (sre(n) ^ 2 + (om - sim(n)) ^ 2)
    Next n
    trans = a0 * g * Sqr(zq / dq)
    Print # 4, Round(f, 6), pp(trans), Round(ph, 6)
  Wend
  Close # 3
  Close # 4
  Print
  Print "          Response was stored in file response.lst"
Else
  Print
  Print "          No file frequencies.lst was found"
EndIf
Print
Input "=============> (9) ok?   ";ok
End
Function atn0(y, x)
  Local phi
  If x = 0
    phi = - pih * Sgn(y)
  Else
'##E#
    phi = Atn(y / x)
  EndIf
  Return phi
EndFunc
Function pp(x)
  ' pretty-print floating-point numbers
  Local ns
  ns = Floor(Log10(x))
  If ns < 5
    Return Str(Round(x, 5 - ns))
  Else
    Return Str(Round((x / 10 ^ ns), 6)) + "E" + Trim(Str(ns))
  EndIf
EndFunc
$DatFile
$ExeName = C:\Dokumente und Einstellungen\Erhard\Desktop\POTSDAM_2013\software\winresp\WINRESP.Exe
$FileDescription = -<({[]})>-
$Comments = Autor: E. Wielandt
$CompanyName = Uni Stuttgart
$LegalCopyright = © 2000 Uni Stuttgart
$LanguageId = 0407
$FileVersion = %d.%d.%d.%d
$ProductVersion = %d.%d.%d.%d
$DatFile