Archive 2024

Section I – R&D strategy forecast

Our Research on 2025 will focus on metering devices tailored for emerging markets. Specifically, our goal is to provide a portable, computer driven phasor measurement unit and power quality analysis solution for Low/Medium voltage (110V-230V to 11kV) at a fraction of the price tag of market leader solutions, which would leverage existing DAC solutions, such as professional USB sound cards instead of costly DAQ platforms or integrated solutions. Most of the R&D effort will be focused towards PC software development. A calibration channel will be provided to take into account hardware induced frequency and phase response variation as well as vendor specific sound cards delays, as well as audio interface (ASIO,Directplay, etc.) delays.

As for financing, we plan to introduce a comprehensive crowdfunding program. We also envision a crowdsourcing initiative for engineers that wish to contribute to this endeavour. This will help in kick starting the project and ensure viable time to market delays. Stay tuned.

A significant part of the project is already underway. We do not start from scratch at the present.

The final implementation will be IEEE Std C37.118.1 compliant.

Section II – 2024 Annual report

2024 Annual report will be made available January 2nd, 2025.

The goal of the following project circuit is to demonstrate the feasibility of audio communication at line levels, with emphasis on voice (male and female) over low voltage (230V line to neutral) power lines.

This project has close to zero usefulness in our current age, besides design practice, for the TX part. It could have some niche use as musical / artistic performance device, or in a recording studio.

We will focus first on the RX device. Since such a device does not inject EMI on power lines, it is not a big no no in terms of use over utility mains. a derivative – and more industrially useful – project could be done on that basis to perform noise and harmonics analysis of the power line environment on the cheap.

Disclaimer

Such a device, specifically the TX part, if constructed, would not be FCC compliant or whatever normative regulation is in place in your country. As such, This article describes a theoretical device and any prototype derived from it should not be coupled to the mains, as it can disrupt the proper function of power meters, inject unwarranted noise on power lines and be hazardous to operate if proper precautions are not taken. Dangerous voltages are present in the device which can cause bodily harm or even death.

Proper handling and safety precautions when designing, prototyping and operating the device are required.

Any practical realization and testing should be done on an isolated power network (on generator power or inverter power). We are not responsible for generator or inverter failure or any other damage or unintended consequences arising from the realization of this project.

Additional disclaimer : EMI mitigation

If coupled to the mains, due to the relative proximity of the voice spectrum to the 50Hz mains frequency, a companion LP line filter device (a ‘reactor’ in electrician speak) at the point where we want to stop audio propagation would be required, and would need large inductors and capacitors. It would, however, add a low impedance path at audio frequencies, and that would absorb part of the energy that is supposed to be received on the RX end, lowering overall SNR.

In the absence of such a filter, the distribution panel and power meter would be exposed to electrical noise in the audio spectrum, which could prevent proper operation of some safety devices.

In case of testing on an isolated network (for instance, on a gasoline or diesel generator), the various filtering stages in the device exhibit a high Q factor and generator frequency stability is paramount in order for the mains hum to be adequately filtered. An inverter generator or an offline inverter test bench is preferable. The following circuit is designed for 50 Hz operation. 60 Hz operation would require offsetting the corner frequency of all filter stages.

Modelling AC mains

We will model an AC mains of 230V Line to Neutral, 50 Hz, and take into account just one phase. A crude model of a local utility 11 kV/ 230V transformer is provided.

Safety

The audio path should be galvanically isolated from mains. For this purpose a 600:600 ohm EI-14 audio transformer is used, and provides adequate bandpass on the audio band. However, such a transformer is limited to 75mW RMS power and has low reactance at 50 Hz. Connecting the primary directly to line and neutral would destroy it immediately. Such a transformer has only 140 Ohms DC resistance on primary and secondary windings and barely much reactance, since its primary inductance (measured at an excitation signal of 1kHz – as an audio standard – is only 230mH). Since magnetizing current is inversely proportional to inductance, it would reach an unsafe level, overheat the transformer windings and melt them.

The solution is to high pass (AC couple) the primary to reduce voltage such as the power dissipated from magnetizing current stays well below 75mW power RMS. This portion of the circuit is fundamentally similar to what is used for PLC (power line communication) : a high pass filter (C3 + R4) plus a pulse transformer. (L1 + L3 // L2 + L4). This forms our first filtering stage.

In our case, instead of a fc at several hundred kHz, it is closer to the power line first odd harmonic (150 Hz) and uses an audio transformer instead of a very low inductance pulse transformer.

Another possible option to reduce the 50 Hz signal at the signal transformer input is to use a capacitive voltage transformer (CVT) which in our case would be a simple capacitive voltage divider. The frequency response is slightly different with a theoretical 6dB resonance boost around 1kHz with the CVT compared to the high pass filter, but with a steeper phase shift due to the resonance, which could give a less natural sound.

Since we are working on audio signals, it could be convenient to make the whole device work on a balanced line. Low cost EI14 600 ohms audio transformers have only two output terminals, so, to emulate a center-tapped transformer, we have to resort to two matched transformers with the secondary wired in series, and we can wire the primary in parallel after the line LP filter. This would get use a 1:2 (pri:sec) voltage ratio and a 40dB/dec rolloff.

The center point of our two secondaries wired in series (L2 + L4) would be our audio ground reference, and also the power ground. it would be bounded to the center tap of a two rail +/-15V power transformer. This other transformer will power op-amps after a rectifier bridge and two 12V regulators (LM7812 + LM7912). The power providing transformer is modeled by (L9 + L7 // L8 + L10). Note that the supplied inductances inductance should be modeled according to the measurement of the device you will use, and that they may skew the operation of the subsequent filters. Which means you will need to adapt the model downstream. In practice however, tuning would be provided by potentiometers to adapt to changing conditions and different upstream component parameters.

Finally the TX device – not modeled – is replaced by a voltage source (V1) stimulated by an audio signal (wave file) containing a voice sample, and with the output characteristics of an audio amplifier : output stage at 8 ohms impedance. At this node the signal is mixed with the mains power (‘line’ node, from the 11kV/230V transformer network)

The output of the first HP (out1/out1b) stage after the audio transformer (L1 + L3 // L2 + L4) exhibits max voltage gain at around 1Khz, at +6dB (due to transformer configuration), a 40dB/dec rolloff high pass, with an attenuation of close to 56dB at 50Hz.

Just downstream of the differential audio line (out1/out1b) transformer output, we add a fully differential buffer op amp so as not to load substantially the audio transformer. a LTC1992 would fit conveniently.

Following that, it would be useful to obtain more 50 Hz rejection by adding a high Q band stop filter tuned at 50 Hz. We will use two filters for each leg of the balanced line, referenced to audio ground. This gets us an additional attenuation of 10.7dB at 50 Hz.

With the following parameters declarations for the filter parameters : (Parametrization permits stepping the filter parameters, to find the best theoretical filter parameters values)

Not all fully differential op amps are capable of performing filtering without ringing artefacts, and that was seen in the simulation. For this purpose, we used two single ended op amps at each stage except for the buffer op amp after the transformer. Calibration should be performed so as to obtain symmetrical and maximum attenuation on both legs, which would require matched op amps and 1% tolerance components, particularly for the band pass and band stop filters.

The other option is to use fixed resistors with the double of the required resistance value, and add a trim potentiometer in parallel that is then fine tuned so that each signal leg performs symmetrically. A // resistor configuration minimises noise during trim pot operation, as wipers may not always contact fully, and particularly during tuning.

For the band stop stage, trim pots would be used across R2,R3,R10 and maybe R11 and the respective resistors on the other leg.

Besides the fully differential buffer, the remaining op amps (for the remaining stages) are LT1037 and performed adequately in the simulation.

Pushing the 50Hz attenuation further.

We will use a complementary circuit that uses a standard 2VA +-15V AC/DC transformer to power up opamps (after a regulator) and to perform an additional function. It can be used, not only to power the device components, but also to provide a 50 Hz signal at a lower voltage than mains, that can be further attenuated with a flat response by a simple voltage divider for it to stay well below the op_amp rail levels.

This signal is not affected by the first pass HP stage. (C3 + R4). It will then be used to obtain better attenuation. However care should be taken to minimise distortion and noise induced from the rectifier diode commutation, by adding a snubber.

We would then perform the opposite of the band stop stage on the first audio path, and add a high Q band pass filter tuned at 50 Hz.

This filtered 50 Hz signal would then be passed to an all pass filter (phase adjust) to make it 180° out-of-phase with our audio signal from the first signal path, and then mixed with a weighted inverting summing amplifier to cancel the 50 Hz signal from our audio line even more. Proper adjustment would use two sets of (coarse/fine) potentiometers, one for the phase and the other for the gain. each potentiomerer could be ganged so as to adjust both legs of the phase adjust and mixer.

A theoretical total attenuation of close to 108dB is obtained at 50 Hz.

Stability of the 50 Hz cancellation.

Since the two mixed signals to achieve active cancellation take two distinct audio paths with different filtering characteristics, the 50 Hz component will be out of phase with respect to each other. The 50 Hz band pass signal path should be brought 180° out of phase before the summing amplifier, and its amplitude should be equal. For this purpose, an all-pass filter is used and tuned to add some phase shift to obtain a 180° phase shift between the two signals.

The following factors have to be taken into account for the cancellation not to drift significantly :

• Frequency stability of the 50 Hz grid : Usually free-running at +- 10mHz around 50 Hz. Above that level, grid primary controls (generators controllers) kick-in to bring back the frequency to 50 Hz, grid-level secondary controls from dispersed PMUs units (synchrophasors) are also used to bring mains frequency to stability, but rare deviations up to 0.1 Hz under heavy power use / power production unbalanced states cannot be ruled out. Since the band stop and band pass filters operate at the corner frequencies, phase shift is in the center of the transition range where d(phi)/df is maximum, and opposite sign between the bandpass and bandstop filters. A slight frequency deviation will make the phase shift drift roughly equal to $$ \phi_{drift} \approx \Delta (f)*(d(\phi_{bandpass}) + d(\phi_{bandstop}))/df $$
• Mains voltage sag and swell : Although both audio paths react to voltage gain changes linearly, if the power transformer used for the 50 Hz bandpass signal path is also used for op amp power supply through a regulator, it can affect loading of the transformer in a non linear manner. This would have a non linear effect on relative amplitude of the signals, and the gain of one leg of the summer amplifier would have to be adjusted. This effect would be more noticable on low power PCB transformers with poor voltage regulation. (where loading creates a higher voltage drop). A high VA rating power transformer with a no load over load voltage ratio close to unity is thus preferable; another option is to use a separate transformer for op-amp power supply. The input impedance of a buck converter can be described by the following formula :
  • 
  • https://e2e.ti.com/support/power-management-group/power-management/f/power-management-forum/787869/snva538-input-impedance-calculation-of-a-buck-converter
  • D is the duty cycle and depends on V(in)/V(out) ratio and is a quite significant term, and quadratic. Frequency effects are negligible for mHz drifts of mains frequency. In our simulation, adding a LM7812 / LM7815 pair of voltage regulators for op-amp supply was sufficient to throw cancellation out of balance.

Sending an audio signal over the power line : Line driver design.

So far, we have only taken into account the receiver circuit design. Now let’s dwelve into driving the power line with an audio signal.

We know that :

• impedance seen at the subscriber level with no load drawing power is quite low, and is mainly dictated by the local 11kV/400V transformer supply energy to the utility customers, as well as any filters located in the power meters or elsewhere that absorb unwanted high frequency signals used in metering or for EMI mitigation. The fc of these filters is mainly dependent on the PLC technology and frequency bands used by power utilities companies.
• Impedance is variable according to power line loading and may include reactive components if inductive loads such as motors are used.
• Broad spectrum noise due to light dimmers, led drivers, switching power supplies and power factor correction circuits is to be expected. This can have an effect on current but on voltage as well.
• Rectifiers induce low frequency harmonics of 50 Hz
• Not all noise sources will spectrally peak in the 20-20kHz audio band.

When using off the shelf line drivers, we have mainly two options : line level driving ( -10dBV consumer or +’dBu professional). +4dBu would be preferable, but the main problem would be output impedance would be still too high to properly drive an AC power line. Line level impedance is typically comprised between 100 to 600 ohms, with the recent trend going to lower impedances. The impedance mismatch is not the only factor to take into account. Line drivers are designed to drive significantly higher load impedances. They are designed to provide a signal, not power. A low impedance load would significantly overstress the output op-amp in terms of current. <current capability of line drivers ?>

The other option being amplifier “hot” signals, such as those designed to drive speakers. Impedance levels are usually 4 to 8 ohms. Galvanic isolation through a high power 1:1 audio signal transformer is highly recomended, plus a high pass filter with the same characteristics as the filter in front of the audio RX transformer, so as not to couple the audio amplifier stage with too high levels of AC 50 Hz.

This would give a first order filter with 20dB/decade slope with a 1: 1 transformer arrangement.

Digital Filter equalization.

Once the RX signal goes through the DAC, in order to obtain a flatter response for audio signals in the audio chain, we have to apply a digital equalization : a low shelf digital filter with a slope of 60dB/decade (we have to take into account the filtering at the TX stage) with a low transition frequency at 20Hz and a high transition frequency at around 730 Hz, where the gain reaches 0dB is necessary. The previously filtered trace 50 Hz signal at the analog level would get a boost, which can be further rejected by a high Q band stop Butterworth digital filter.

Going full duplex.

It would be nice to use a single device for both RX and TX, and filter the TX signal generated at the same endpoint, leaving only the reciprocal signal heard from the remote endpoint. This is known as sidetone reduction. On our RX side a 180° out of phase copy of the signal sent to the TX driver would be combined in the RX audio path. Depending on the placement, different equalization curves would be required on the TX signal copy to account for the return signal filtered characteristics. Gain should be adjusted in tandem (at the TX send level and sidetone mixing level) through either carefully calibrated ganged potentiometers or through an AGC setup. Detecting the absence of sidetone through analog means is not an easy task. digital correlation filters on the other hand could be used.

When interfacing a Geiger-Muller device with a MCU or a Raspberry pi, there are mainly two options available :

• Use of a high level device, which transfers counter information digitally such as cps, cpm, tube type, stored memory or whatever information the device holds. Some devices may also be able to inform the host of a pulse detection using a low level signaling mode, and perform as a low level device. High level devices may be more or less sophisticated and allow for calibration, unit conversions, dosimetry, or changing settings depending on the tube model and their efficiency characteristics. High count compensation due to dead time may or may not be done, but any serious device should perform such processing. The compensation formulas are dependent on the designer knowledge of the matter and their proprietary choices.

• Use of a low level device. Such a device is an analog front end that supplies high voltage to the GM tube, and outputs detected events as pulses on a GPIO pin. The short pulse is usually stretched to a fixed time pulse, with sharp rising and falling edges. It should essentially output a square wave signal, that is easily processed by subsequent A/D stages that are not part of the device itself. Signal output is commonly done via GPIO header pins and/or 3.5 audio jacks. Most also feature a buzzer for direct audio pulse output.

This article is focused on the low level devices and their interfacing with a digital counting platform such as the Raspberry Pi, as it offers maximum processing flexibility.

Note that for high CPM processing, digital latency effects starts to be significative when using a non real-time OS platform due to operating system overhead. Part of this article will deal with mitigation of these unwanted effects.

GM Tube basics.

There are plenty of resources explaining the physics of the GM Tube. We won’t cover them in detail except for the important issue of dead time and how this skew the measured count rate at high levels of ionizing radiation (high count rates). This issue is critical as the device underestimates the count at these regimes. That is, the device is not linear.

The GM Tube is supplied a high enough voltage to work in the Geiger plateau, such as an ionization event triggers a Townsend avalanche and a significant charge is displaced between anode and cathode, and this is seen by the amplifier (either a transistor or an op-amp) as a voltage pulse. This charge displacement is equivalent to a lowering of the impedance of the GM Tube, which induces a voltage drop from the nominal high voltage, as the high voltage source is far from ideal, it can’t source much current and has added resistors. The GM tube and resistors on the high voltage path perform as a voltage divider (with a parallel path to ground, the GM tube – with its impedance lowering during the Townsend avalanche process)

This voltage drop is conditioned to a low voltage (with reference to ground), buffered, and the pulse is stretched to several hundreds of microseconds to a couple of milliseconds, depending on the design and/or tuning.
Pulse stretching usually involves stages such as Schmitt triggers and 555 timers

Dead time from the analog front end.

Usually the analog front end circuitry processing exhibits a non-paralyzable dead-time characteristic.

dead time : it means that after a pulse is detected by this front end, any subsequent pulse is not detected during the interval spanning from the registered pulse to the expiration of the dead-time (the pulse stretched “HIGH’ state time, set usually by an analog timer such as a 555.)

non-paralyzable : it means that any townsend avalanche event occuring during that dead-time interval. does not re-trigger the dead time, or at least, it should not. If it does, it would be a poor design choice as dead-time from any processing stage should be minimized as much as possible, and a paralyzable dead-time may extend this duration.

GM Tube intrinsic dead time.

As for the GM tube Townsend avalanche process, it has an intrinsic dead time characteristic that is short but still significative. After a sucessful ionizing event detection, one that triggers a full Townsend avalanche, the tube is unable to detect any other event during the dead time. This dead time is usually in the order of 80 to 200 µs, and is mostly paralyzable : Any event ocurring during that dead time, resets that intrinsic dead time to some extent. For sake of simplicity, we will treat it as a fully paralyzable dead time. In reality, most of the literature agrees that this intrinsic dead time falls in between a paralyzable and non paralyzable model.

We make the supposition that the dead time is a function of the time separation of two subsequent (or more) ionization events, with the limit reaching the nominal dead time as event separation in the time domain goes toward infinity.
The precise model of a Townsend avalanche and its practical implementation in a commercial GM tube is outside of the scope of this article.

Thus it follows that the dead-time nominal values from the GM Tube and those of the analog pulse stretch front end are important parameters for count rate compensation.

Analog frontend deadtime measurement is quite straightforward in a low count (low background radiation) environment, and can be done by plugging the GM meter jack output to a sound card set at the highest sampling capability, ideally 96ks/s or more. Or better, by using an oscilloscope.

In our case, the dead time was measured at around 2.2ms average on this commercial device.

As for the intrinsic dead time of the GM Tube, some are reported in the datasheets. Determination of that dead time is not an easy task, in part because the ionization events exhibit a probabilistic distribution (modeled as a Poisson process) – That is, it is hard to control the precise moment a ionization event will trigger a Townsend avalanche.

Furthermore, GM Tubes have an efficiency parameter: Of all the gamma flux that intercepts the section area of the tube (usually the tube length times the tube diameter, and taking into account the mass absorption characteristics of the tube wall due to limb effects if one wants a really precise characterization), not all gamma flux will interact with the tube gas electrons, the main factors to take into account are the finite – and small – tube volume, low gas density, and type of gas used as it influences ionization behaviour. Note that GM tubes are sensitive to x-ray (with quite low efficiency) gamma rays, and charged particles such as alpha, beta+ and beta-. Non charged particles such as neutrons are detected indirectly, by a nuclear process that give rise to charged particles when a neutron is captured or absorbed by the GM tube gas. Neutron sensitive GM tubes may use BF3 gas (Boron Tri-Fluoride)

GM Tube efficiency

Common beta/gamma small glass tubes exhibit an efficiency of around 0.03

Moreover, this efficiency varies with the gamma ray energy. It is well known that GM Tubes under-perform in the X-ray spectrum. Flattening of the energy/efficiency response curve can be done by using filters (such as thin metal sheets) that partially occult the tube, mainly for gamma counting.

As for tubes that exhibit a better sensitivity in the alpha and beta detection mode, they usually are made with walls that have a lower extinction coefficient for alpha or beta particle, such as mica.

Remember also that GM tubes are not proportional counters, scintillators, or particles chambers. They do not convey information about the source of the ionization event energy, as the Townsend avalanche “saturates” the detector.

True count rate estimation based on GM tube dead time and analog front end dead time

Based on the research of Muller in “Dead Time problems” – Nuclear instruments and methods 112th issue (1973), (1) The model we will use is stated on page 56 (d) equation (56), as it accounts for two dead times in series, one paralyzable from the GM tube, and one non-paralyzable from the analog front end (pulse stretcher)

we can see that the analog dead time due to pulse shaping into a square wave is the limiting factor as for high counting rates. Thus, it is preferable to tune, in the design phase, the 555 timing to achieve the shortest pulse length that does not trigger spurious counts by the A/D stage. A/D stages typically register the rising edge of the pulse. As the amplitude of the pulse is not of any particular interest, a high frequency, low bit resolution A/D design is preferable, such as a sigma/delta ADC, (ADC considerations go beyond our scope, as we are limited by the BCM GPIO hardware on the Raspberry Pi) and good noise immunity practices, such as using shielding for the whole assembly and a shielded signal cable grounded at one end, or twisted pair, while maintaining the shortest possible signal path from the GM front end to the ADC. This is to ensure good signal integrity and low stray loop inductance. A very short pulse could induce ringing and be registered as multiple events when the level lingers between the low/high logic boundary. However, there are diminishing returns for high count precision when the analog front end dead time approaches the intrinsic GM tube dead time, that is in the order or 80 to 200 µs. These time constants are well within the 555 timer capabilities, and proper design should take care of ringing artifacts. As for the hardware limitations of Raspberry pi pulse sensing, they are in the tens of ns range. Which means that lowering analog front-end dead time for ADC sensing to around 100 µs should not be a hardware limiting factor. As for the software GPIO / ISR library, it is recommended to use one that minimizes delay. Some libraries are implemented as linux daemons that may allow better management of high count rates. Best is to compare performance with a signal generator between the standard Python Raspberry Pi GPIO libraries that perform better with high pulse frequency (in the order of 10 kHz).

Python code processing could be the bottleneck, if care is not taken to benchmark deque() performance at 10 kHz pulse rate, particularly if the deque() stores the pulse timing with calls to Python time libraries.

Some GPIO pulse processing libraries store pulse timestamp information, in this case this overhead could be discarded.

Also, registering the falling edge may be useful to assess proper operation of the GM analog front end. Absence of a falling edge in a timely manner before the next rising edge could signify that the signal is stuck in a high state, and should display a malfunction and/or high count warning.

It is also preferable to use a micro controller or full fledged miniature computer such as a Raspberry pi with adequate processing power, which translates into a fast CPU clock and more than one core (for computers) to decrease ISR (interrupt service request) burden to a minimum in high count environment.

In the case of a Raspberry Pi, and when using Python, the following guidelines should be followed :

• Minimum amount of code in the ISR.
• Ideally it should use a deque() for pulse registering, simply appending the pulse to the deque in the ISR, and exiting the ISR.
• Assessing the need for specialized GPIO libraries for pulse counting, such as those that involve a daemon, when pulse timestamping with the highest precision is a design requirement. In that case however, inter-process communication (IPC) delay is a factor to take into account for overall responsiveness.
• A separate thread on another core should perform post processing such as logging into the file-system or a database, or count rate compensation calculations.
• It is preferable to use a low level language such as C++ instead of Python, and precise benchmarking should be done using a function generator generating pulses with a comparable duty time to the GM pulse shaper backend, with increasing CPM frequency, to assess the digital latency induced and ISR responsiveness in high CPM environments. This way, the influence of the A/D backend and code performance can be precisely factored in for final count up-rating & device calibration.

Pulse time-stamping

In our design the deque() stores a timestamp for each pulse. Pulse time-stamping is a niche requirement and may be of use to triangulate the source of a sudden burst of gamma energy, although air attenuation coefficient factor coupled with poor efficiency of GM tubes, would render such an endeavour tricky. As for alpha and beta particle radiation, their detection on fixed position sensors for environmental monitoring in weather stations is dependent on wind and atmospheric currents drift that carry contaminated dust and fallout, which operate on timescales large enough to not require precise time stamping.

For EAS (extensive air shower) research arising from high energy cosmic particles, as well as study of dark lightning generated terrestiral gamma-ray flashes (TGF) a scintillation detector would be a better choice, as they have a better quantum yield.

A GPS module is a good investment in a project of this kind as it allows not only Geo tagging of events, but also precise timestamping due to inherent time synchronisation features of GPS.

Alternatively, low quality synchronisation of nodes that perform event timestamping may use NTP or higher precision NTP protocols. In any case, precise node time synchronisation using NTP requires symmetric network packet processing (no asymmetric routing – this creates different propagation delays upstream and downstream). These propagation delay uncertainties increase substantially when the time source is several router hops away, and also depend on the network traffic load induced delay.

If time synchronisation is performed through air (such as using Lora) all radio induced propagation delays and bitrate induced delays have to be factored in.

Count up-rating using two dead times in series model with time constants t1 and t2.

In our project, we will use Muller’s derivation (1) p56. (d)

$$ R = \frac{\rho}{(1-\alpha)x} + e^{\alpha x} $$

$$ x = \rho t_{2} $$

$$ \alpha = \frac{t_{1}}{t_{2}} $$

$$ \rho $$

being the true count rate estimation, and

$$ R $$

being the measured count

$$ t_{1} $$

being the (paralyzable) dead time of the GM Tube, and

$$ t_{2} $$

being the (non-paralyzable) dead time of the analog frontend pulse shaper.

Let’s introduce the other well known models accounting for a single dead time system:

the non paralyzable model :

$$ R = \frac{\rho}{1 + t\rho} $$
the paralyzable model :

$$ R = \rho e^{(-\rho t)} $$

Since the unknown is the corrected count (rho), we need to use the inverse function of these models, regardless of the model, compound, paralyzable or non paralyzable.

The paralyzable function inverse expression requires the use of the W0 and W1 Lambert function, Math helpers in Python such as scipy allow straightforward calculation of the Lambert W0 and W1 branches, albeit with some computational burden.

The compound t1 and t2 in series requires numerical methods such as the secant method. Which would only increase the computational burden.

In the case of a Raspberry Pi, since RAM and storage are not an issue, and the problem is not multivariate since t1 and t2 are constants. We advise to compute the functions models, and use a reverse table lookup for fast determination of the corrected count. Scipy propose linear and higher order interpolation mechanisms, which would have a lower computational burden than root finding.

Calibration and experimental confirmation of the count correction model.

Disclaimer : For this part, access to a medium activity point source > 100 kBq, with low uncertainty is preferable. A low activity source would not push the GM counter at count rates where there is significant deviation from the true rate, and thus not give enough data to properly test the models. Depending on your location, point sources are exempt of declaration at different ceilings. To complicate things further, several regulations may overlap such as national / European Union, and providers of sources may or may not be available in your country. a source in the 100 kBq pretty much requires a license anywhere in the world. One option may be to use the source of a third party at the site of the third party licensed metrology / research lab. As always, exposure to a medium activity source may be harmful depending on exposure time and proper handling and shielding precautions are required.

Besides the point source, a lead plated rectangular cross section channel spanning from the point source to the GM Tube is preferable, such a device function is not to serve as a collimator, but rather to prevent reflections of gamma rays leading to constructive / destructive interference, as some papers suggest. It seems quite remarkable however that such effects interfere significantly with the measurement, given the very low refractive index of most materials in the gamma spectrum. They would however shield the detector somewhat from background radiation and other sources that may be found in a lab at relative proximity.

Note that a Cs137 nucleide decays either to a stable 56Ba137 nucleide through Beta- decay, with a probability of 5.4%, or to an excited (56Ba137m) state, also through Beta- decay, with a probability of 94.6%. When reverting to the ground state 56Ba137, some of the 56Ba137m nucleides emit a 0.6617 MeV gamma photon. Of all Cs137 decays that form the measure of the total activity of the sample supplied by the manufacturer, 85.1 % yield gamma photons. This has to be factored in, besides manufacturer supplied source percent uncertainity, and exponential law derating of the activity of the sample, using the supplied date of manufacture, and the Cs137 half-life of 30.05 years.

Beta- radiation has a high linear attenuation coefficient in air, and would skew the measurements when the jig is very close to the tube through their contribution in the final measure. A standard food grade aluminium foil is sufficient to filter them out to a negligible level, while leaving almost all gamma ray photons unscathed. Nevertheless we have factored in attenuation from air and aluminium through air and aluminium for both Beta- and gamma based on litterature and NIST database data.

The source is then placed on a jig powered by a linear actuator, so that it moves freely inside the channel along the x axis. the center of the point source disc should be on the x axis and intercept the center of mass of the GM tube. The x axis should be perpendicular to the point source disc and GM tube.

The raspberry Pi platforms registers pulse and operates the linear actuator. A model with position feedback is required for accurate jig position determination, after manual position calibration, or alternatively, a time of flight (TOF) / LIDAR sensor module should be used. This solution should also be subjected to calibration beforehand.

The GM calibration helper for a Cs-137 point source Python script is in alpha stage of development but the bulk of the work is done.

The workflow of the script starts with the following parameter inputs :

• Point source activity, supplied in either Bq or Ci units.
• The activity_unit enumeration is used to set up the unit of activity
• The radioisotope used is hard-coded to Cs137, and it’s decay scheme is factored in to take into account only Gamma activity.
• Decay % that give rise to gamma photons
• Cs137 half-life as a constant
• Source date of manufacture
• Aluminium mass extinction coefficients for Beta and gamma radiation (using Cs137 decay energies)
• Air mass extinction coefficients for Beta and gamma radiation (using Cs137 decay energies)
• Distance range (min,max) of the source relative to the GM tube center of mass on the jig y-axis
• GM Tube geometry : length and diameter.
• GM tube intrinsic dead time estimation (or from datasheet)
• GM tube efficiency estimation (or from datasheet)

On that basis, we will first calculate the mean path length from the point source to the tube as a function of y axis distance. For this we will use basic trigonometry and integrate over the tube length, all the ray paths from the source. This mean path length will be used to factor in the linear attenuation coefficient of air for gamma radiation at 0.611 MeV energy. Beta radiation is assumed to be filtered up to an insignificant amount, but we can still calculate the attenuation based on aluminium mass extinction coefficient and foil thickness to check for filtering efficiency of beta radiation.

The compound attenuation of air plus aluminium filter is calculated for the gamma radiation.

The radiant flux from the point source received by the GM tube is then calculated. The standard law

$$ Flux =\frac{P}{4 \pi r^{2}} $$

assumes a spherical irradiated surface. since the cross section of the GM Tube is a plane, and taking into account that $$ GMtubewidth \ll GMtubelength $$, we will perform a single integration along the tube cylinder axis to get the corrected number of photons $$ R_{net} $$ crossing the tube section, instead of a double integration.

$$ R_{totvac} \simeq 2\times\int_{0}^{\frac{GMtubelength}{2}}\frac{P_{net}}{4 \pi (x^2 + r^2)} cos(arctan(\frac{x}{r}))dx \hspace{2mm} \times GMtubewidth $$

Which simplifies to :

$$ R_{totvac} \simeq 2\times\frac{P_{net}\frac{GMtubelength}{2}}{4\pi r \sqrt{(\frac{GMtubelength}{2})^{2} + r^{2}}} \times GMtubewidth $$

Also,

$$ P_{net} $$ is the gamma source activity. It has to take into account the percent of decays giving rise to gamma photons, and the reduction of activity since the source date of manufacture.

$$ P_{net} = 0.851 P_{nom} e^{-\frac{ln(2)}{t_{1/2}} (t_{1} – t_{0})} $$

$$ t_{1} $$ is the time of the calibration

$$ t_{0} $$ is the time of source manufacture.

and, accounting for attenuation :

$$ R_{tot} \simeq 2\times att_{\gamma}\frac{P_{net}\frac{GMtubelength}{2}}{4\pi r \sqrt{(\frac{GMtubelength}{2})^{2} + r^{2}}} \times GMtubewidth $$

To calculate the air attenuation of the gamma photons, we will use an approximation that has an expression in terms of elementary functions.

A gamma photon path length ‘l’ from the point source to any point on the x axis (axis of the GM tube represented as a cylinder) can be approximated as the length of the hypotenuse (neglecting paths that strike the cylinder off axis) :

$$ l = \sqrt{x^{2} + r^{2}} $$

The average path length is :

$$ l_{avg} = {\frac{2}{GM tube length}} \int_{0}^{\frac{GM tube length}{2}} \sqrt{x^2 + r^2}dx $$

Which yields the following expression :

$$ l_{avg} = \frac{{r^{2} \operatorname{arsinh}\left(\frac{GMtubelength}{2r}\right)}}{GMtubelength} + \frac{\sqrt{({\frac{GMtubelength}{2}})^{2} + r^{2}}}{2} $$

The attenuation factor applied to the flux :

$$ att_{\gamma}\simeq e^{-\mu_{air} l_{avg} -\mu_{Al}w} $$

$$ w $$ being the Aluminium foil thickness

$$ \mu_{air} $$ being the linear attenuation coefficient of air at 0.661 MeV gamma energy.

$$ \mu_{Al} $$ being the linear attenuation coefficient of Aluminium at 0.661 MeV gamma energy.

A small fraction of these photons will trigger an ionization of the GM_tube gas medium. These will form the measured tube event count per second. $$ R_{net} $$ At low count rates, dead time effects are negligible, and provided sufficiently long measurement times, the tube efficiency alpha can be determined :

$$ \alpha = \frac{R_{net}} {R_{tot}} $$

Resources

Note : For up to date code, please check https://www.github.com/rodv92/GMobile

This is the calibration helper Python script draft.

import scipy as sp
import math as m
import numpy as np
import time as time
from datetime import datetime
import matplotlib.pyplot as plt
from enum import Enum


# using S.I units (unless specified otherwise in comment immediately following 
# declaration)

class activity_unit(Enum):
    BQ = 1 #Bq
    CI = 2 #Ci

activity_unit = Enum('activity_unit', ['BQ','CI'])



# assuming a Cs-137 source, and a collimator made of thick lead of the same width as the GM tube and the same height as the tube, with the point source at the e# ntrance of the channel


P_nom = 370e-3 # source activity (unit determined by used_activity_unit)
used_activity_unit = activity_unit.BQ
# it's better to use a high activity source in order to drive the count at very high levels for dead-time effects to be come noticeable
# the lower the analog frontend dead-time that results in reliable A/D conversion, the higher the activity of the source is required.


gamma_yield = 0.851 # ratio of decays that give rise to gamma photons

source_date_of_manufacture = "2022-12-01T19:00:00Z"
dt_source_dom = datetime.fromisoformat(source_date_of_manufacture)
dt_now = datetime.now()
source_age = (dt_now - dt_source_dom).total_seconds()
half_life = 30.05 # Cs137 half-life in years
half_life = half_life*365*24*60*60 # Cs137 half-life converted to seconds 

µ1 = 0.000103 # linear attenuation coefficient of air for the gamma energy of 667 kEv - emmited by Ba137m
#µ2 = 0.1 # linear attenuation coefficient of beta shield for the gamma energy of 667 kEv. We assume that the beta shield filters beta particles to a negligible amount. 
Al_filter_thickness = 0.001 # cm  0.001 cm = 10µm standard food grade aluminium foil thickness


max_distance = 0.3
distance = 0.15 # distance from point source to center of detector. (meters)
min_distance = 0.05
# GM tube axis is normal to line from point source to GM tube center.

GM_tube_length = 0.075 # glass envelope length only (meters)
GM_tube_diameter = 0.01 # glass envelope diameter (meters)
GM_tube_alpha = 0.0319 # tube efficiency - not all incoming photons trigger an ionization event.
alpha = GM_tube_alpha

GM_tube_detection_cross_section_area = GM_tube_length*GM_tube_diameter
GMT_dcsa = GM_tube_detection_cross_section_area # variable alias

GM_tube_dead_time = 80e-6 # tube intrinsic (charge drift time) dead time
GMT_det = GM_tube_dead_time

if (used_activity_unit == activity_unit.BQ):
    pass
elif (used_activity_unit == activity_unit.CI):
    P_nom *= 3.7e10 # convert to Bq


#derate activity based on % decay gamma yield and source age.
P_net = gamma_yield*P_nom*exp(-(m.log(2)/half_life)*source_age)

#calculate mean path length of gamma photons reaching the tube.
mean_path = (1/2)*(distance**2 + (GM_tube_length/2)**2)**(1/2) + (distance**2)*m.arcsinh*(GM_tube_length/(2*distance))/GM_tube_length
#calculate attenuation from linear attenuation coefficient
gamma_att_air = m.exp(-µ1*mean_path)


Al_density = 2.7 # g/cm3
#Al_mass_att_beta_cs137 = 15.1 #cm2/mg
Al_mass_att_beta_cs137 = 15.1e3 #cm2/g beta+/- attenuation for Cs137 emitter https://doi.org/10.1016/j.anucene.2013.07.023
# not used in subsequent calibration formulas, assuming that beta is filtered to an insignificant amount. standard food grade aluminium foil of 10µm thickness gives an attenuation factor in the order of 1e-18

µ2 = Al_mass_att_beta_cs137*Al_density #cm-1
beta_att_Al = m.exp(-Al_filter_thickness*µ2)
# To check the effectiveness of beta filtering.
# not used in subsequent calibration formulas, assuming that beta is filtered to an insignificant amount. standard food grade aluminium foil of 10µm thickness gives an attenuation factor in the order of 1e-18

Al_mass_att_gamma_Ba137m = 7.484e-2 #cm2/g gamma attenuation for Aluminium at 667 kEv
µ3 = Al_mass_att_gamma_Ba137m*Al_density

gamma_att_Al = m.exp(-Al_filter_thickness*µ3)

#calculate total gamma attenuation from air and Al filter.

gamma_att_total = gamma_att_air*gamma_att_Al




import GMobile as GM


# CALIBRATION step 0

# Estimates the analog front end dead-time by timing the pulse width while the GM Tube is exposed to background radiation

background_measure_pulsewidth_total_pulses = 60 # time to spend in seconds measuring pulse width
background_measure_pulsewidth_max_fails = 20 # time to spend in seconds measuring pulse width
rise_timeout_ms = 20000
fall_timeout_ms = 20


(pulsewidth,stdev) = GM.measurePulseWidth(background_measure_pulsewidth_total_pulses,background_measure_pulsewidth_max_fails,rise_timeout_ms,fall_timeout_ms)


# CALIBRATION step 1

#Measure the background CPM with the collimating assembly, but the source removed and far away.
# call main.py and average CPM over specified background_acquire_time in seconds


background_acquire_time = 600 # time to spend in seconds acquiring background radiation levels after first 60 sec of acquisition.


chars = None
chars = input("Step 1 - acquiring background radiation cpm during " + background_acquire_time + " seconds. Please put the source as far away as possible. press ENTER to start")
while chars is None:
    chars = input("Step 1 - acquiring background radiation cpm during " + background_acquire_time + " seconds. Please put the source as far away as possible. press ENTER to start")
    

GM.SetupGPIOEventDetect() # sets up the GPIO event callback

s = 0
cpm_sum = 0
while (s < background_acquire_time):
    cpm = GM.process_events(False,False) # This call should take exactly one second.
    if (cpm != -1):
        cpm_sum += cpm
        s += 1

cpm_background = cpm_sum/background_acquire_time

def model1_estimated_GM_CPM(true_count,t1,t2): # Muller's serial t1 paralyzable dead time followed by t2 non paralyzable dead time model
    alpha = t1/t2
    x = true_count*t2
    corrected_cpm_1 = true_count/((1-alpha)*x + m.exp(alpha*x))
    return corrected_cpm_1

def model2_estimated_GM_CPM(true_count,t1):
    corrected_cpm_2 = true_count*m.exp(-true_count*t1)
    return corrected_cpm_2

def model3_estimated_GM_CPM(true_count,t2):
    corrected_cpm_3 = true_count/(1+true_count*t2)
    return corrected_cpm_3


def movejig(position):

    #TODO : linear actuator positioning code
    err = 0
    
    return err #  err = 0 : actuation OK

def efficiency_step(distance=0.15,movestep=0.005,efficiency_placing_time=600,efficiency_stab_time=60,last_secs_stab_time=60,min_cpm_efficiency_cal=8*cpm_background,max_cpm_efficiency_cal=16*cpm_background):

    # CALIBRATION step 2

    print("Step 2 : automated GM tube efficiency calculation")
        
    s = 0
    s_stab = 0
    cpm_stab = []
    while (s < efficiency_placing_time):
        #print("Step 2.1 - efficiency computation: you have " + efficiency_placing_time + "seconds to put the source at a distance to obtain a reading between " + min_cpm_efficiency_cal + " and " + max_cpm_efficiency_cal + " cpm. press ENTER when in range")
        #print("the timer will start counting down after the first 60 seconds have elapsed - it will then exit the step if cpm is stabilized in range for a whole " + efficiency_stab_time + " secs, and get the cpm average for the last " + last_secs_stab_time + " secs. Do not move the source during that time")
        cpm = GM.process_events(False,False) # This call should take exactly one second.
        if (cpm != -1): # first 60 seconds have elapsed.
            cpm_sum += cpm
            s += 1
            if(cpm > min_cpm_efficiency_cal and cpm < max_cpm_efficiency_cal): # in range.
                s_stab += 1
                cpm_stab.append(cpm)
            elif(cpm > max_cpm_efficiency_cal): # out of range high. reset stabilized countrate time counter
                distance += movestep
                movejig(distance)
                s_stab = 0
                cpm_stab = []
            elif(cpm < min_cpm_efficiency_cal): # out of range high. reset stabilized countrate time counter
                distance -= movestep
                movejig(distance)
                s_stab = 0
                cpm_stab = []
                
            
            print("cpm:\t" + cpm)
            print("stabilized_time:\t" + s_stab)
            print("distance:\t" + distance)
            
            
            if s_stab >= efficiency_stab_time:
                #distance = input("Please input distance in meters from GM_tube to source at stabilized reading")
                cpm_stab = cpm_stab[-last_secs_stab_time:]
                cpm_stab_avg = sum(cpm_stab)/last_secs_stab_time
                break

    cpm_efficiency_calc = cpm_stab_avg - cpm_background

    #calculate mean path length of gamma photons reaching the tube.
    mean_path = (1/2)*(distance**2 + (GM_tube_length/2)**2)**(1/2) + (distance**2)*m.arcsinh*(GM_tube_length/(2*distance))/GM_tube_length
    #calculate attenuation from linear attenuation coefficient
    gamma_att_air = m.exp(-µ1*mean_path)

    #calculate total gamma attenuation from air and Al filter.
    gamma_att_total = gamma_att_air*gamma_att_Al


    flux = gamma_att_total*2*GM_tube_width*(P_net*(GM_tube_length/2))/(4*m.pi*distance*((GM_tubelength/2)**2 + distance**2)**(1/2))
# gamma flux in photons.s^-1 crossing the GM tube.  Accounting for planar cross section of GM Tube (instead of solid angle) and attenuation from air and beta filter of gamma photons.
    theoretical_cpm = 60*flux # assuming GM tube efficiency of 1, All photons would give rise to ionization events inside the GM tube

    if(s != efficiency_placing_time):
        efficiency = cpm_efficiency_calc/theoretical_cpm
        print("computed efficiency:\t" + efficiency)
        return((efficiency,distance))
    else:
        print("efficiency calibration failed.")
        return((-1,distance))

# Compute efficiency of detection at a count rate sufficiently high enough above background but not as high as dead time effects become significant.
efficiency_placing_time = 600
efficiency_stab_time = 180
last_secs_stab_time = 60
min_cpm_efficiency_cal = 8*cpm_background
max_cpm_efficiency_cal = 16*cpm_background


chars = None
chars = input("Step 2.0 - efficiency computation: please put the source at " + distance + "cm from the GM tube, normal to the tube, withint the collimator. press ENTER to start")

movejig(distance) # initial position
(efficiency,distance) = efficiency_step(distance) # gets efficiency, -1 if failed, and jig to source distance

while (efficiency == -1):
    print("efficiency calculation step failed. repeating")
    (efficiency,distance) = efficiency_step(distance)


#STEP 3 : record countrate while stepping the source jig towards the GM tube
step3_wait_time = 120
step3_last_seconds_measure = 60
step3_last_seconds_measure = min(120,step3_last_seconds_measure) # ensure the total seconds we sample is lower than step3_wait_time
x_distance = np.arange(max_distance,min_distance,-0.005)

y_cpm_m = np.empty(len(x_distance)) # array of average of count rate for each measurement sampling
y_cpm_m[:] = np.nan

y_std_m = np.empty(len(x_distance)) # array of standard deviation of count rate for each measurement sampling
y_std_m[:] = np.nan

y_cpm_t = np.empty(len(x_distance)) # array of theoretical cpms derated with GM tube efficiency
y_cpm_t[:] = np.nan

y_cpm_tm = np.empty(len(x_distance)) # array of theoretical cpms derated with GM tube efficiency and dead time effects
y_cpm_tm[:] = np.nan


distance_cpm_avg_m = np.stack([x_distance,y_cpm_m],axis=1) # tabular data for cpm (average) as a function of source/GM tube distance
distance_cpm_std_m =  np.stack([x_distance,y_std_m],axis=1) # tabular data for cpm (std dev) as a function of source/GM tube distance

distance_cpm_t = np.stack([x_distance,y_cpm_t],axis=1) # tabular data for cpm (theoretical), derated by GM tube efficiency as a function of distance
distance_cpm_tm = np.stack([x_distance,y_cpm_tm],axis=1) # tabular data for cpm (theoretical), derated by GM tube efficiency and accounting for dead time effects


distance = max_distance # retract actuator to minimum to get longest source to GM tube distance.
if not (movejig(distance)): # no actuation error
    idx = 0
    while(distance > min_distance):
        if(movejig(distance)):
            break # actuation error, break loop
        distance -= 0.05
        cpm_stab = [] 
        #TODO : reset deque() in GMobile after jig move, to get rid of the inertia induced by the sliding window sampling
        while(s < step3_wait_time):
            if(s >= (step3_wait_time - step3_last_seconds_measure)):
                cpm_stab.append(GM.process_events(False,False)) # This call should take exactly one second.
            else:
                time.sleep(1)
        cpm_avg = np.average(cpm_stab)
        cpm_std = np.std(cpm_stab)
        distance_cpm_avg_m[idx][1] = cpm_avg
        distance_cpm_std_m[idx][1] = cpm_std
        
        
        #calculate mean path length of gamma photons reaching the tube.
        mean_path = (1/2)*(distance**2 + (GM_tube_length/2)**2)**(1/2) + (distance**2)*m.arcsinh*(GM_tube_length/(2*distance))/GM_tube_length
        #calculate attenuation from linear attenuation coefficient
        gamma_att_air = m.exp(-µ1*mean_path)
        #calculate total gamma attenuation from air and Al filter.
        gamma_att_total = gamma_att_air*gamma_att_Al

        flux = gamma_att_total*2*GM_tube_width*(P_net*(GM_tube_length/2))/(4*m.pi*distance*((GM_tubelength/2)**2 + distance**2)**(1/2))

        # gamma flux in photons.s^-1 crossing the GM tube. Accounting for planar cross section of GM Tube (instead of solid angle) and attenuation from air and beta filter of gamma photons.
        theoretical_cpm_eff = 60*flux*efficiency # theoretical cpm (derated with GM tube efficiency estimated in step 2.0)

        distance_cpm_t[idx][1] = theoretical_cpm_eff
        distance_cpm_tm[idx][1] = model1_estimated_GM_CPM(theoretical_cpm_eff) # theoretical cpm from above with dead time compensation




def compare_cpm_measured_theoretical(theoretical,measured):

    devsum = 0
    devratiosum = 0

    if (len(theoretical) != len(measured)):
        return (-1,-1)
    
    for idx in range(0,len(theoretical)):

        devratiosum += abs((measured[idx][1] - theoretical[idx][1])/measured[idx][1])
        devsum = (measured[idx][1] - theoretical[idx][1])**2

    mape = devratiosum/len(theoretical)
    return (devsum,mape)

def plotcurves(curve_x, curve_y1, curve_y2, color_curve_1, color_curve_2, label_curve_1, label_curve_2):
    

    plt.plot(curve_x, curve_y1, color_curve_1, label=label_curve_1)
    plt.plot(curve_x, curve_y2, color_curve_2, label=label_curve_2)

    plt.xlabel('distance (m)')
    plt.ylabel('count rate (cpm)')
    plt.legend()
    plt.grid()
    plt.show()

print(compare_cpm_measured_theoretical(distance_cpm_tm,distance_cpm_avg_m))
plotcurves(distance_cpm_t[:,0],distance_cpm_avg_m[:,1],distance_cpm_tm[:,1],'r-', 'g-', "measured","theoretical_compensated")

This is the main GM counter interfacing script, with Mariadb logging, individual pulse timestamping Geotagging, and Lora telemetry (work in progress)

# coding: utf-8


import RPi.GPIO as GPIO
import scipy as sp
import signal
import sys
import time
import datetime
from collections import deque 
import struct
import statistics

# Module Imports
import mariadb
#from mariadb.connector.aio import connect
import sys
import re

# Lora Imports

from SX127x.LoRa import *
#from SX127x.LoRaArgumentParser import LoRaArgumentParser
from SX127x.board_config import BOARD

def signal_handler(sig, frame):
    GPIO.cleanup()
    sys.exit(0)


def pulse_detected_callback(channel):
    global pulse_events
    pulse_events.append(time.time() - start_time_epoch)


def dms2dec(dms_str):
    
    
    dms_str = re.sub(r'\s', '', dms_str)
    
    sign = -1 if re.search('[swSW]', dms_str) else 1
    
    numbers = [*filter(len, re.split('\D+', dms_str, maxsplit=4))]

    degree = numbers[0]
    minute = numbers[1] if len(numbers) >= 2 else '0'
    second = numbers[2] if len(numbers) >= 3 else '0'
    frac_seconds = numbers[3] if len(numbers) >= 4 else '0'
    
    second += "." + frac_seconds
    return sign * (int(degree) + float(minute) / 60 + float(second) / 3600)

# Connect to MariaDB Platform
try:
    conn = mariadb.connect(
    #conn = connect(
        user="gmobile_user",
        password="gmobile_passwd",
        host="127.0.0.1",
        port=3306,
        database="gmobile"

    )
except mariadb.Error as e:
    print(f"Error connecting to MariaDB Platform: {e}")
    sys.exit(1)

# Get Cursor
cur = conn.cursor()

start_time_epoch = time.time()
pulse_events = deque()
cps = 0
cpm = 0
#process_pulses = True

#paralyzable model
# m = n*exp(-nt)
# using lambert(w) 
#y =  w * exp(w)
#-nt = w

#y = -nt * exp(-nt)

#-nt = lambert(y)
#n = -lambert(y)/t

#scipy.special.lambertw(z, k=0, tol=1e-8)


#while True:
#    GPIO.wait_for_edge(PULSE_GPIO, GPIO.FALLING)
#    print("Button pressed!")
#    sp.special.lambertw(z, k=0, tol=1e-8)

PULSE_GPIO = 16
lat = dms2dec("51°24'04.30\"N")
long = dms2dec("30°02'50.70\"E")


def setupGPIOEventDetect():


    GPIO.setmode(GPIO.BOARD)
    GPIO.setup(PULSE_GPIO, GPIO.IN)
    GPIO.add_event_detect(PULSE_GPIO, GPIO.RISING, callback=pulse_detected_callback)


def removeGPIOEventDetect():

    GPIO.remove_event_detect(PULSE_GPIO)


def measurePulseWidth(total_pulses,max_fails=20,rise_timeout_ms=20000,fall_timeout_ms=20):

    removeGPIOEventDetect()
    s = 0
    fail = 0
    pulsewidths = []

    while(s < total_pulses and fail < max_fails):
        channel = GPIO.wait_for_edge(PULSE_GPIO,GPIO.RISING,timeout=rise_timeout_ms) # assuming there is at least one pulse registered every rise_timeout_ms
        if (channel is None):
            fail += 1
            continue
        pulsestart = time.time_ns()
        channel = GPIO.wait_for_edge(PULSE_GPIO,GPIO.FALLING,timeout=fall_timeout_ms) # assuming the pulse width is less fall_timeout_ms
        if (channel is None):
            fail += 1
            continue
        pulsewidths.append((time.time_ns() - pulsestart)/1000) # save pulsewidth value in µs
        s+= 1 # sucessful pulsewidth measure

        #this pulse timing method may be problematic if the real time separation between RISING and FALLING edge is shorter than execution time between the two
        #wait_for_edge calls. thread should be set to high priority. Overall, dead time will be overestimated.

    setupGPIOEventDetect()
    return (pulsewidths[int(len(pulsewidths)/2)],statistics.pstdev(pulsewidths)) # returns mean pulsewidth value in µs and standard deviation.



# LORA INIT

#BOARD.setup()
#BOARD.reset()
##parser = LoRaArgumentParser("Lora tester")







"""
class mylora(LoRa):
    
    global current_data
    
    def __init__(self, verbose=False):
        super(mylora, self).__init__(verbose)
        self.set_mode(MODE.SLEEP)
        self.set_dio_mapping([0] * 6)

    def on_rx_done(self):
        BOARD.led_on()
        #print("\nRxDone")
        self.clear_irq_flags(RxDone=1)
        payload = self.read_payload(nocheck=True )# Receive INF
        print ("Receive: ")
        mens=bytes(payload).decode("utf-8",'ignore')
        mens=mens[2:-1] #to discard \x00\x00 and \x00 at the end
        print(mens)
        BOARD.led_off()
        if mens=="INF":
            print("Received data request INF")
            time.sleep(2)
            print ("Send mens: DATA RASPBERRY PI")
            self.write_payload([255, 255, 0, 0, 68, 65, 84, 65, 32, 82, 65, 83, 80, 66, 69, 82, 82, 89, 32, 80, 73, 0]) # Send DATA RASPBERRY PI
            self.set_mode(MODE.TX)
        time.sleep(2)
        self.reset_ptr_rx()
        self.set_mode(MODE.RXCONT)

    def on_tx_done(self):
        print("\nTxDone")
        print(self.get_irq_flags())

    def on_cad_done(self):
        print("\non_CadDone")
        print(self.get_irq_flags())

    def on_rx_timeout(self):
        print("\non_RxTimeout")
        print(self.get_irq_flags())

    def on_valid_header(self):
        print("\non_ValidHeader")
        print(self.get_irq_flags())

    def on_payload_crc_error(self):
        print("\non_PayloadCrcError")
        print(self.get_irq_flags())

    def on_fhss_change_channel(self):
        print("\non_FhssChangeChannel")
        print(self.get_irq_flags())

    def start(self):          
        while True:
            self.reset_ptr_rx()
            self.set_mode(MODE.RXCONT) # Receiver mode
            while True:
                pass;

    def send(data):
        data = [255, 255, 0, 0] + data + [0] 
        self.write_payload([data]) # Send DATA
        time.sleep(2)
        self.reset_ptr_rx()
        self.set_mode(MODE.RXCONT)

lora = mylora(verbose=False)
#args = parser.parse_args(lora) # configs in LoRaArgumentParser.py

#     Slow+long range  Bw = 125 kHz, Cr = 4/8, Sf = 4096chips/symbol, CRC on. 13 dBm
lora.set_pa_config(pa_select=1, max_power=21, output_power=15)
lora.set_bw(BW.BW125)
lora.set_coding_rate(CODING_RATE.CR4_8)
lora.set_spreading_factor(12)
lora.set_rx_crc(True)
#lora.set_lna_gain(GAIN.G1)
#lora.set_implicit_header_mode(False)
lora.set_low_data_rate_optim(True)

#  Medium Range  Defaults after init are 434.0MHz, Bw = 125 kHz, Cr = 4/5, Sf = 128chips/symbol, CRC on 13 dBm
#lora.set_pa_config(pa_select=1)



assert(lora.get_agc_auto_on() == 1)

try:
    print("START")
    lora.start()
except KeyboardInterrupt:
    sys.stdout.flush()
    print("Exit")
    sys.stderr.write("KeyboardInterrupt\n")
finally:
    sys.stdout.flush()
    print("Exit")
    lora.set_mode(MODE.SLEEP)
BOARD.teardown()
"""

   
signal.signal(signal.SIGINT, signal_handler)

def process_events(log=False,lorasend=False):
    global cps
    global cpm
    nowtime = time.time()
    prune_before_time = nowtime - start_time_epoch - 60.0
    for pulse in list(pulse_events):
        if(pulse < prune_before_time):
            pulse_events.popleft()
    cps = len(pulse_events)/60.0
    cpm = len(pulse_events)
    print(cpm)
    
    if not(int(time.time()) % 60) and prune_before_time > 0:
        # log last minute cpm to db. wait at least 60 sec from start
        # to get steady state data.
        print("last minute cpm:" + str(cpm))
        if(lorasend):        
            epoch_ms = int(time.time()*1000.0)
            buffer_data = struct.pack("<Q", epoch_ms) # pack epoch_ms into byte array little endian
            buffer_data += struct.pack("<L", cpm) # append cpm into byte array little endian
            #lora.send(buffer_data)
            print(buffer_data)
        if(log):
            sql = "INSERT INTO data_cpm (timestamp_utc,count_per_minute,coordinates) VALUES(%s, %s, ST_PointFromWKB(ST_AsBinary(POINT(" + str(lat) + "," + str(long) + ")), 4326))";
            val = (str(datetime.datetime.now()),str(cpm))
            cur.execute(sql, val)
            conn.commit()
    
    processing_delay = time.time() - nowtime
    print(processing_delay)
    proc_delay_mul = int(processing_delay)
    time.sleep(1 + proc_delay_mul - processing_delay)
    # ensure isochronous sampling, add n skip seconds in case of the block taking a really long time...
    # asynchronous mariadb cursor should prevent this occurence.
    if(prune_before_time > 0): # 60 sec have elapsed
        return cpm
    else:
        return -1

    
#while process_pulses:
#    process_events()
#signal.pause()