Kurt's Weblog

I went through the cracks data a little more this morning and verified that all the layback points I used have there speed from within a few shots. And thankfully that is the case. Here is the results table with all 52 entries.

#   1     2     3     4       5      6        7
# shot wireout vert layback angle1 speed  shot_offset_for_speed



525106 264.6 79.8168 252.2745299 17.56 3.8 2
528154 240.8 85.0668 225.273788 20.69 3.8 17
595405 218.5 62.5668 209.3505327 16.64 5.0 0
597545 254.7 82.8168 240.8598506 18.97 3.8 0
598722 254.7 79.0668 242.1167717 18.09 4.2 0
605053 254.7 76.0668 243.0759798 17.38 3.7 17
633320 185.6 73.0668 170.6124343 23.18 4.0 0
670983 258.6 76.0668 247.1594666 17.11 3.9 0
683526 278.3 89.5668 263.4932226 18.77 4.0 0
689024 256.6 73.8168 245.7532096 16.72 4.1 0
711636 191.8 82.3668 173.2135972 25.43 4.2 0
713188 244.8 82.3668 230.5271139 19.66 4.3 0
741727 212.4 62.5668 202.9757511 17.13 4.8 0
745538 272.5 62.5668 265.2199946 13.27 4.5 0
760000 272.5 73.8168 262.3115134 15.72 3.8 0
815797 185.6 73.0668 170.6124343 23.18 3.8 0
821962 240.8 93.3168 221.9833661 22.80 3.6 0
844793 216.5 73.0668 203.797676 19.72 3.7 0
850749 270.5 76.0668 259.5844601 16.33 4.2 0
856363 305.2 76.0668 295.5687432 14.43 4.0 0
913976 214.4 64.0668 204.6040203 17.39 4.2 0
918322 258.6 79.0668 246.2161675 17.80 4.1 0
921739 295.6 91.0668 281.2226839 17.94 4.3 0
946083 185.6 72.3168 170.9316835 22.93 4.1 0
952200 314.8 84.3168 303.2980667 15.54 4.2 0
1289000 185.6 73.8168 170.2892834 23.44 3.9 0
1306641 242.8 79.0668 229.565418 19.00 4.0 0
1355961 210.3 65.5668 199.8176287 18.17 4.2 0
1360480 242.8 85.8168 227.1284149 20.70 4.0 0
1386006 191.8 76.8168 175.7453249 23.61 3.8 0
1387816 248.7 85.8168 233.4248634 20.19 4.1 0
1396015 248.7 72.0168 238.0446818 16.83 3.8 0
1405500 204.1 63.0168 194.1280323 17.98 4.0 0
1433000 195.6 64.0668 184.8101868 19.12 4.9 0
1433557 262.6 76.8168 251.1133992 17.01 4.5 0
28067 191.8 72.3168 177.6443651 22.15 3.6 0
29794 278.3 92.5668 262.4543342 19.43 3.8 0
37030 278.3 96.3168 261.101444 20.25 3.6 0
39700 305.2 91.8168 291.0613599 17.51 3.7 0
44800 272.5 72.3168 262.7290057 15.39 3.9 0
60575 128.4 56.5668 115.2681966 26.14 3.7 0
61542 310.9 84.3168 299.2482034 15.74 3.6 0
223345 177.2 61.0668 166.3450809 20.16 3.6 0
229052 195.6 79.8168 178.5739019 24.08 3.9 0
263470 170.9 67.0668 157.1905033 23.11 3.7 0
267038 266.5 82.0668 253.5493844 17.94 3.6 0
309101 305.2 85.8168 292.8865255 16.33 3.3 0
313592 258.6 83.5668 244.7254583 18.85 4.8 0
316027 208.2 74.5668 194.3888689 20.99 3.7 0
383600 258.6 76.0668 247.1594666 17.11 3.6 0
394351 220.6 61.0668 211.9792583 16.07 3.7 0
418500 256.6 80.4918 243.6485792 18.28 4.2 0

After verifying the data, I created a little least square fit program to use the results table.

#!/usr/bin/env python
########################################
# FIT a line to the data
dataFile = open('results.dat','r')
dataFile.readline()
dataFile.readline()
dataFile.readline()
wireoutAngle=[]
speedAngle=[]
for line in dataFile.readlines():
    fields = line.split()
    wireoutAngle.append((float(fields[2-1]),float(fields[5-1])))
    speedAngle.append  ((float(fields[6-1]),float(fields[5-1])))
#
from Scientific.Functions.LeastSquares import leastSquaresFit
def linear(parameters,values):
    a,b = parameters
    x = values
    return (a*x + b)
#
initial = (-.1,32)
param,error = leastSquaresFit(linear,initial,wireoutAngle)
print "  ----- WIREOUT TO ANGLE ----- "
print param,error
print 130, 130 * param[0] + param[1]
print 300, 300 * param[0] + param[1]
#
initial = (0,22)
param,error = leastSquaresFit(linear,initial,speedAngle)
print "  ----- SPEED TO ANGLE ----- "
print param,error
print 3.4, 3.4 * param[0] + param[1]
print 4.8, 4.8 * param[0] + param[1]

That made a simple linear fit to the data. The results of a run look like this:

./cracks_cruise_layback_fit.py
  ----- WIREOUT TO ANGLE ----- 
[-0.051842754250418684, 31.465098411287943] 194.768779108
130 24.7255403587
300 15.9122721362
  ----- SPEED TO ANGLE ----- 
[-1.703222454552509, 25.827653874339632] 415.848028016
3.4 20.0366975289
4.8 17.6521860925

After creating two data files for gnuplot, I end up with plots that have the best fit line.





The code for the processing of the excel data from Jenna is in cracks_cruise_layback.py.html (raw python: cracks_cruise_layback.py)

Finally a quick calculation from the Santa Barbara Cruise...

ipython
In [1]: from math import *
In [2]: from units import rad2deg



In [3]: rad2deg(asin(30.75/50))
Out[3]: 37.951920289246438  # Wire angle seen



In [4]: 50 * (-0.051842754250418684) + 31.465098411287943
Out[4]: 28.87296069876701   # Wire angle predicted

I went through the cracks data a little more this morning and verified that all the layback points I used have there speed from within a few shots. And thankfully that is the case. Here is the results table with all 52 entries.

#   1     2     3     4       5      6        7
# shot wireout vert layback angle1 speed  shot_offset_for_speed



525106 264.6 79.8168 252.2745299 17.56 3.8 2
528154 240.8 85.0668 225.273788 20.69 3.8 17
595405 218.5 62.5668 209.3505327 16.64 5.0 0
597545 254.7 82.8168 240.8598506 18.97 3.8 0
598722 254.7 79.0668 242.1167717 18.09 4.2 0
605053 254.7 76.0668 243.0759798 17.38 3.7 17
633320 185.6 73.0668 170.6124343 23.18 4.0 0
670983 258.6 76.0668 247.1594666 17.11 3.9 0
683526 278.3 89.5668 263.4932226 18.77 4.0 0
689024 256.6 73.8168 245.7532096 16.72 4.1 0
711636 191.8 82.3668 173.2135972 25.43 4.2 0
713188 244.8 82.3668 230.5271139 19.66 4.3 0
741727 212.4 62.5668 202.9757511 17.13 4.8 0
745538 272.5 62.5668 265.2199946 13.27 4.5 0
760000 272.5 73.8168 262.3115134 15.72 3.8 0
815797 185.6 73.0668 170.6124343 23.18 3.8 0
821962 240.8 93.3168 221.9833661 22.80 3.6 0
844793 216.5 73.0668 203.797676 19.72 3.7 0
850749 270.5 76.0668 259.5844601 16.33 4.2 0
856363 305.2 76.0668 295.5687432 14.43 4.0 0
913976 214.4 64.0668 204.6040203 17.39 4.2 0
918322 258.6 79.0668 246.2161675 17.80 4.1 0
921739 295.6 91.0668 281.2226839 17.94 4.3 0
946083 185.6 72.3168 170.9316835 22.93 4.1 0
952200 314.8 84.3168 303.2980667 15.54 4.2 0
1289000 185.6 73.8168 170.2892834 23.44 3.9 0
1306641 242.8 79.0668 229.565418 19.00 4.0 0
1355961 210.3 65.5668 199.8176287 18.17 4.2 0
1360480 242.8 85.8168 227.1284149 20.70 4.0 0
1386006 191.8 76.8168 175.7453249 23.61 3.8 0
1387816 248.7 85.8168 233.4248634 20.19 4.1 0
1396015 248.7 72.0168 238.0446818 16.83 3.8 0
1405500 204.1 63.0168 194.1280323 17.98 4.0 0
1433000 195.6 64.0668 184.8101868 19.12 4.9 0
1433557 262.6 76.8168 251.1133992 17.01 4.5 0
28067 191.8 72.3168 177.6443651 22.15 3.6 0
29794 278.3 92.5668 262.4543342 19.43 3.8 0
37030 278.3 96.3168 261.101444 20.25 3.6 0
39700 305.2 91.8168 291.0613599 17.51 3.7 0
44800 272.5 72.3168 262.7290057 15.39 3.9 0
60575 128.4 56.5668 115.2681966 26.14 3.7 0
61542 310.9 84.3168 299.2482034 15.74 3.6 0
223345 177.2 61.0668 166.3450809 20.16 3.6 0
229052 195.6 79.8168 178.5739019 24.08 3.9 0
263470 170.9 67.0668 157.1905033 23.11 3.7 0
267038 266.5 82.0668 253.5493844 17.94 3.6 0
309101 305.2 85.8168 292.8865255 16.33 3.3 0
313592 258.6 83.5668 244.7254583 18.85 4.8 0
316027 208.2 74.5668 194.3888689 20.99 3.7 0
383600 258.6 76.0668 247.1594666 17.11 3.6 0
394351 220.6 61.0668 211.9792583 16.07 3.7 0
418500 256.6 80.4918 243.6485792 18.28 4.2 0

After verifying the data, I created a little least square fit program to use the results table.

#!/usr/bin/env python
########################################
# FIT a line to the data
dataFile = open('results.dat','r')
dataFile.readline()
dataFile.readline()
dataFile.readline()
wireoutAngle=[]
speedAngle=[]
for line in dataFile.readlines():
    fields = line.split()
    wireoutAngle.append((float(fields[2-1]),float(fields[5-1])))
    speedAngle.append  ((float(fields[6-1]),float(fields[5-1])))
#
from Scientific.Functions.LeastSquares import leastSquaresFit
def linear(parameters,values):
    a,b = parameters
    x = values
    return (a*x + b)
#
initial = (-.1,32)
param,error = leastSquaresFit(linear,initial,wireoutAngle)
print "  ----- WIREOUT TO ANGLE ----- "
print param,error
print 130, 130 * param[0] + param[1]
print 300, 300 * param[0] + param[1]
#
initial = (0,22)
param,error = leastSquaresFit(linear,initial,speedAngle)
print "  ----- SPEED TO ANGLE ----- "
print param,error
print 3.4, 3.4 * param[0] + param[1]
print 4.8, 4.8 * param[0] + param[1]

That made a simple linear fit to the data. The results of a run look like this:

./cracks_cruise_layback_fit.py
  ----- WIREOUT TO ANGLE ----- 
[-0.051842754250418684, 31.465098411287943] 194.768779108
130 24.7255403587
300 15.9122721362
  ----- SPEED TO ANGLE ----- 
[-1.703222454552509, 25.827653874339632] 415.848028016
3.4 20.0366975289
4.8 17.6521860925

After creating two data files for gnuplot, I end up with plots that have the best fit line.





The code for the processing of the excel data from Jenna is in cracks_cruise_layback.py.html (raw python: cracks_cruise_layback.py)

Finally a quick calculation from the Santa Barbara Cruise...

ipython
In [1]: from math import *
In [2]: from units import rad2deg



In [3]: rad2deg(asin(30.75/50))
Out[3]: 37.951920289246438  # Wire angle seen



In [4]: 50 * (-0.051842754250418684) + 31.465098411287943
Out[4]: 28.87296069876701   # Wire angle predicted

Jenna has done just over 50 layback calculations for one cracks cruise off of the eastern seaboard. I just did a short analysis of the relationship between wireout, layback and speed. Here are the resulting graphs. These results may possibly have a couple bad points that I still need to get rid of.





The impression I get is that the wire angle is more dependent on the amount of wire in the water than the speed of the ship.

The data file that results from my processing looks like this.

#   1     2     3     4       5      6
# shot wireout vert layback angle1 speed



525106 264.6 79.8168 252.2745299 17.56 3.8
528154 240.8 85.0668 225.273788 20.69 3.8
595405 218.5 62.5668 209.3505327 16.64 5.0
597545 254.7 82.8168 240.8598506 18.97 3.8
598722 254.7 79.0668 242.1167717 18.09 4.2
605053 254.7 76.0668 243.0759798 17.38 3.7
633320 185.6 73.0668 170.6124343 23.18 4.0
670983 258.6 76.0668 247.1594666 17.11 3.9
683526 278.3 89.5668 263.4932226 18.77 4.0
689024 256.6 73.8168 245.7532096 16.72 4.1
711636 191.8 82.3668 173.2135972 25.43 4.2
713188 244.8 82.3668 230.5271139 19.66 4.3
741727 212.4 62.5668 202.9757511 17.13 4.8
745538 272.5 62.5668 265.2199946 13.27 4.5
760000 272.5 73.8168 262.3115134 15.72 3.8

I need to cleanup the code some before I share it, but it is getting there. Right now, my worry is that I had to remove any point from the shot-speed table where I had two shotpoint numbers that were the same. The current version of the code will keep scanning until it find the next shot with speed and use that speed to match up with the wireout and layback.

I was thinking about using the python random class to make some simple figures that need random particles dispersed in a uniform manner. Since random needs a seed, I figured a good way to give it one is to use /dev/random. It is easy to use the python struct class to turn the random data that /dev/random generates into values of a certain type. Here the 'i' tells struct that the random data is a 4 byte integer.

#!/usr/bin/env python



import struct
class DevRandom:
    """http://www.python.org/doc/current/lib/module-struct.html"""
    def __init__(self):
        self.devrand = open ('/dev/random','r')
        return
    def getInt4(self):
        string = self.devrand.read(4)
        intval = struct.unpack('i',string)
        return intval[0]



r = DevRandom()
for i in range(100000):
    print r.getInt4()

Just released pmag-kds-py-0.7. This adds the xdr.py from my previous post.

• xdr.py.1.html - documents the command line program
• xdr.html - documents the code

Thanks to Wei Wei for explaining how to get through the software.

xrd.py now lives in pmag-kds-py. It reads the data from the XRD machine in the SIO Analytical Facilities.

xrd.py --cps-file=4-63b-cps.txt --show --no-cps --esd

I have just released pmag-kds-py 0.6. There will be no more pmag-py releases from me. pmag-kds-py contains only my python code. It is up to Lisa as to if, when, and how any of pmag-kds-py gets merged into pmag-py and when pmag-py will be released.

Just to reiterate, pmag-kds-py is a use at your own risk piece of software without support of any kind. If you are going to use it in your research, I expect you to verify that the calculations in here are correct. I am sure there are errors in code.

http://schwehr.org/software/pmag-kds-py/

schwehr.org/software/pmag is no more.

This is pretty exciting (other than the not available for Mac OSX and not an open source db)... Oracle to offer free database

  The database heavyweight on Tuesday is expected to announce the beta
  release of Oracle 10g Express Edition (Oracle Database XE), which
  will be generally available by the end of the year. It is targeted
  at students, small organizations and software vendors that could
  embed the Oracle database with an application. 



  The latest edition is the same as other databases in Oracle's lineup
  but is limited in usage. It can only run servers with one processor,
  with 4GB of disk memory and 1GB of memory. 

I still don't know why oracle would be any better than postgres, but it is still great.

Question: "What do you do for your job?"

Answer: "I shoot laser beams at mud."

Today Jason and I did our first X-ray Diffraction (XRD) measurements today. I know this is not a huge deal for all of you who do XRD all the time, but it is new for us.

SIO's Analytical Facility has a Scintag XDS 2000, which looks like this:



Bruce Deck setup the instrument for us and ran the software. Here are the settings that we used on the front.



I created the quick look sample by scraping the dried cube with a razor blade to fill up a 1 cm square patch around the dimple in the glass slide. The slides are to small for the sample holder, so it has to go in diagonally.



Once the system is closed up with the lead glass shield down and the interlock disabled, it is time to run the sample. The measurements are done in degrees which are actually the emitter and receiver moving up in angle above the sample. Here is the emitter:



The software to control the system is DMSNT on windows... it unfortunately requires a floppy disk to copy the data off. Here is the software running.





Once the software is done running the system through your range (for me 2 to 50°) it can help you try to identify the materials present. Here is the original and baseline corrected data with an illite standard below. The middle graph are the peaks that the machine has found. The machine currently has a 0.2° offset.

Here is pretty much the simplest case of a pdflib_py wrapper that I am making. The API that comes with pdflib for python feels like it came straight out of a C programmers mouth. This is still pretty basic. It only writes some text. It takes an awful lot of work just to get that done in the original API. I am a little more able to cope with a simple object oriented API.

#!/usr/bin/env python
import pdflib_py
import sys
class PDF:
    """OO Wrapper for pdflib_py"""
    def __init__(self,filename):#,size='letter'):
        
        self.pdf = pdflib_py.PDF_new()
        if -1 == pdflib_py.PDF_open_file(self.pdf,filename):
            print "Error: " + pdflib_py.PDF_get_errmsg(self.pdf) + "\n"
            sys.exit(ams.EXIT_FAILURE)
        pdflib_py.PDF_begin_page(self.pdf,595, 842)        # FIX: make page setup sane
        # Make a default font
        self.font = pdflib_py.PDF_load_font(self.pdf,"Helvetica-Bold","iso8859-1","")
        pdflib_py.PDF_setfont(self.pdf,self.font, 12)
    def set_info(self,creator,author=None,title=None):
        if None!=creator: pdflib_py.PDF_set_info(self.pdf, 'Creator', creator)
        if None!= author: pdflib_py.PDF_set_info(self.pdf, 'Author',  author)
        if None!=  title: pdflib_py.PDF_set_info(self.pdf, 'Title',   title)
    def show(self,string,pos=None):
        if None!=pos: pdflib_py.PDF_set_text_pos(self.pdf,pos[0],pos[1])
        pdflib_py.PDF_show(self.pdf,str(string))
    def __del__(self):
        import pdflib_py # bad ref counting in pdflib_py?
        pdflib_py.PDF_end_page(self.pdf)
        pdflib_py.PDF_close(self.pdf)
        pdflib_py.PDF_delete(self.pdf)
        
import sys, datetime
if __name__ == '__main__':
    pdf = PDF("foo.pdf")
    pdf.set_info(sys.argv[0],'Kurt Schwehr','clay flocs')
    pdf.show(datetime.date.today(),(25,25))

Wait! Why am I trying to mess with pdflib? If I use matplotlib or gnuplot, I can still write in python and I do not need to deal with the pdflib_py API. Feel free to take this idea and run with it, but I'm not going to right now. I was just thinking about making a pdf figure and my brain went back to Feb of this year.

How does python's destruction system work? I am cookin up a little pdflib wrapper and am having some strange behavior on cleanup. Yes, I do know that pdflib is up to 6.0.2 and fink only has 5.0.3. Soon, I promiss.

It seems like the pdf object was already garbage collected when __del__ was being reached, so the pdf cleanup is freaking out. This little python snippet proves taht objects are not eaten by the great GC monster before my __del__ code gets run.

#!/usr/bin/env python
class destroy:
    def __init__(self):
        self.a = 1
    def __del__(self):
        print self.a
        del(self.a)
        print self.a
# Try it...
d = destroy()
del(d)

I expect the 2nd print of self.a to throw an exception. And that is what happens.

./destroy.py 
1
Exception exceptions.AttributeError: "destroy instance has no 
  attribute 'a'" in <bound method destroy.__del__ of 
  <__main__.destroy instance at 0x38ce90>> ignored

Aha!!! I figured the bugger out. It was not my pdf object that was gone when my code got to __del__, rather pdflib_py was already gone. Doing another import did the trick, but is that okay? It feels REALLY WRONG to do this sort of thing. Shouldn't the pdf object be holding up the reference count on the import pdflib_py instance? Perhaps the pdflib python interface is missing a reference count? Or is this how python is supposed to work?

    def __del__(self):
        print 'running __del__', type(self.pdf), self.pdf
        import pdflib_py
        pdflib_py.PDF_end_page(self.pdf)
        pdflib_py.PDF_close(self.pdf)
        pdflib_py.PDF_delete(self.pdf)

I am leaning towards a missing refcount in the pdflib_py library. I tried this little snippet and it prints out 3.1415... just fine.

#!/usr/bin/env python
import math
class destroy2:
    def __init__(self):
        return
    def __del__(self):
        print math.pi
d = destroy2()

BTW, No luck trying to figure out what is wrong with mapserver. It seems to get confused and think that my mac needs -arch i386 and -arch ppc. Not gonna work like that. How to make the linker confused.

From this month's OGC News

...
  GlobeXplorer has over 600Tb of imagery content and they add as much
  as 1Tb more each week. This is mostly commercial imagery data
  (including over 2M square kilometers of Digital Globe content). It
  also includes other remotely sensed data sources and ancillary map
  data overlays that can be overlaid onto images (e.g. flood zones and
  data parcel data.) Most of the imagery consists of 6-inch and
  sub-meter resolution data from around the world, and includes some
  data at 3-inch resolution. GlobeXplorer has a very popular portal
  for individual consumers who might want to buy, for example, a high
  resolution image of their neighborhood, but most of their images are
  served to business customers via Web services. They serve over 2
  million maps/day.
...

I just noticed that mapserver is in fink. When did that happen? Oh, June of this year. Wahoo! Now we need gdal python support and OpenEV. Accept mapserver is not quite building on 10.3.9. Whoops.

OpenIOOS - Oceanographic data interchange.

I am hoping to do some X-Ray Diffraction work soon and have two sources of info that look to be good starting points. Thanks to Jill and Mike.

X-Ray Diffraction and the Identification and Analysis of Clay Minerals 
by 
Duane M. Moore
Robert C. Reynolds, Jr
1989

and

Data Report: Normalization Factors for Semiquantitative X-Ray
Diffraction Analysis, with Application to DSDP Site 297, Shikoku
Basin
by
Michael B. Underwood,2 Nandini Basu,2 Joan Steurer,2 and Swati Udas

I am sitting here typing in the data from last weeks SEM run.

QUESTION: How can a sample have -1.24% carbon? Ahem.

This is a note to anyone designing SEM software. Please, oh please put these features into your software. The more info the data files, the better. This stuff could go in the tiff tags, in an XML file, or flat file even. Then I can write a little pdflib program to assemble a composite image. With some EDX machine readable data with position and time, I can even lay in the spot on the pdf.

• Stage position including height
• The scale of the image or width
• Focusing info
• Any other tunable beam parameters
• SEM mode (environmental, low vac, or high vac)
• time and date of the acquisition system
• Oh, and enable lossless tiff compression!
• Don't even let people write jpeg's. What is the point of saving bad data?
• Build a secure system that can publish the data to just me over the net. The SIO WinXP boxes have been force complete off the network. What a waist of CD-R's.
• Why do I have to explicitly store images. I would like a mode that grabs frames at a slow frequency, then let me explicitly mark the ones I need to pull out... Then we could do a pan/scan of most of a sample and build a large Fledermaus or iv model of the surface.

For the EDX mode, exporting a report to a word document is very cute, but a royal pain. A machine readable format is critical.

• A simple ASCII table of the spectra would rock
• Put in the time and stage position for linking to the SEM images
• Include the peak locations for the elements
• Include both the found and missing elements that the detector can sense (knowing that there is no detectable S is important too)
• Why not make EDX acquisition be automatic. The moment the SEM stops moving, start sampling. Then let me know how many samples at the current spot.

Be warned that these are notes from my first use an SEM. Maybe some of these features really are there and I just don't know about them.

We also need a free and open source system for best guess mineral identification. I bet a little WEKA setup could do a good matching job.

Back when I was working on MER, I had the IVS (Fledermaus folks) make a custom setup for microscope work. Both sides need to converge before the average user can just make a 3D model from an SEM and do analysis right in that tool. This is the kind of stuff that I wanted (and still want) to do with a Viz type architecture. Let me build a Rose diagram and/or 3D probability density map in the 3D environment quickly. Tag locations, annotate them, and then drop in additional data. Finally, let me make a publication quality figure and HD quality movie. All without killing myself in the process.

Oh the joy of magnetic units... pmag 0.5 is out:

• Makefile: build docs correctly
• agfm.py: Add getSlope(data), --save-cgs
• agfm.py: High field susceptebility should be correct now
• ams.py: Hext now can normalize on sample volume and is in SI
• units.py: AGFM slope to SI vol normalize susceptibility. Ouch
• units.py: Oe2T() and T2Oe()

http://schwehr.org/software/pmag/

Up until now, the Hext class in ams.py has just returned the values from Lisa's k15_hext. That means KLY-2 raw data which is in microSI will end up generating a bulk susceptibility also in microSI. I am switching the sense of the Hext class to convert to microSI unless you tell it not to. I have only used the KLY-2 at SIO, so I have no idea what other Kappabridges return. Here is a quick example of the new values in my Santa Barbara database. The value for this sample used to be 145.95. The database desperately needs a table that tracks units, but I won't get to do that anytime soon.

sqlite bpsio04.db \
    'select ams.bulksusc,ams_geo.bulksusc \
        from ams,ams_geo \
        where ams.samplename="bp04-1gw-s2-026" 
          and ams.samplename=ams_geo.samplename' | tr '|' ' '



0.000145950012 0.000145950012

This change will be in the pmag-py 0.5 release.

This is fantastic news!!!! Doxygen now supports python. Thanks to Dennis Benzinger for pointing this out to me.

Doxygen features now lists python as the second bullet.

Dennis pointed me to the doxygen formulas page. The best part is that it uses LaTeX (well that is great for Lisa and I), so I don't have to learn MathML right now.

Directions on doxygen for python

There is also pyfilt that mushes python code with doxygen tags into a C file that doxygen can use for documentation.

Jenna is working up some more layback info for the chirp. The hope is to have a lookup table of ship speed to a range of layback angles. Here is a little bit on the commands I am using to generate the velocity input table for the work. Hopefully, the results of this will end up in layback.py

The first task is to remove any bad xstar segy files. Neal's chirp always writes an ASCII header at the front that starts with a company entry like this:

C1  COMPANY:

A quick grep will list those files without the company entry and move them out of the way. They will confuse segysql.py and cause troubles.

mkdir Bad
for file in `grep --files-without-match 'C1  COMPANY' *.xstar`: do
  mv $file Bad
done

This gives me all the munch files:

ls Bad/



line183.sgy     line350.sgy     line396.sgy     line400.sgy     line406.sgy
line336.sgy     line351.sgy     line397.sgy     line403.sgy     line407.sgy
line347.sgy     line352.sgy     line398.sgy     line404.sgy     line408.sgy
line349.sgy     line353.sgy     line399.sgy     line405.sgy     line409.sgy

The most likely cause of the errors are bad tape reads or attempting to read past the end of data on a tape. But don't quote me on that since I am not an expert with tape drives!

Next, we need a short script to generate databases from the xstar files. Since sqlite 2 gets really slow with monster databases, I am going to break them into 5 groups of 100 files. Here is make-db.bash:

#!/bin/bash



 segysql.py --noisy -d xstar -f lines001-099.db line[0-9].sgy line[0-9][0-9].sgy
 segysql.py --noisy -d xstar -f lines100-199.db line1[0-9][0-9].sgy
 segysql.py --noisy -d xstar -f lines200-299.db line2[0-9][0-9].sgy 
 segysql.py --noisy -d xstar -f lines300-399.db line3[0-9][0-9].sgy 
 segysql.py --noisy -d xstar -f lines400-499.db line4[0-9][0-9].sgy 

Once the database is created, here is a sample line used to make the lookup table.

sqlite lines001-099.db 'select shotpoint,speed/10.0,x/3600.0,y/3600.0 \
                        from segytrace' \
                        | tr '|' ' ' > 001-099.shot-speed

Which generates stuff that looks like this with speed in knots.

bzcat 001-099.shot-speed.bz2 | head
4729 3.6 -74.6377777777778 36.8233333333333
4729 3.6 -74.6380555555556 36.8233333333333
4730 3.6 -74.6380555555556 36.8233333333333
4730 3.6 -74.6380555555556 36.8233333333333
4731 3.6 -74.6380555555556 36.8233333333333
4731 3.6 -74.6380555555556 36.8233333333333
4732 3.6 -74.6380555555556 36.8233333333333
4732 3.6 -74.6380555555556 36.8233333333333
4733 3.6 -74.6380555555556 36.8233333333333
4733 3.6 -74.6380555555556 36.8233333333333

Which looks like it is off the eastern seaboard of the US like it is supposed to be.

I really wish that I had a button to enable and disable Javascript and animated gifs on a per page basis. Those annoying animated adds really crank up the CPU useage and eat through the battery on my laptop.

A little more work with Illustrator provided an easy solve for the maximum zoom with these super detailed images - make the lowest resolution image much larger! All I have to do is scale up the first image by 400% and start from there. Now I have no problem with 10000x images. I also overlayed one of the comparison graphs from the EDX and highlighted where the measurements were taken and the spot size.



Wacky Towns: Paducah, Possum Trot and Monkey's Eyebrow. I've seen the exit for Boring, OR. which is not too far from Drain, OR.

Coastal Bluffs Provide More Sand To California Beaches Than Previously Believed

...
  Ashford said decades-old photographs of the Southern California
  coast taken from the ground and the air also have documented the
  steady pace of erosion. However, he said the photographs lack the
  precision and accuracy of the laser scanning technique called LIDAR,
  an acronym for light detection and ranging. Ashford said the 3-D
  maps generated by LIDAR permitted Young to calculate the
  unexpectedly high volume of bluff material that has fallen onto
  beaches during the study period.
...
  At the wave washed western edge of the campus, Neal Driscoll, a
  geology professor at UCSD's Scripps Institution of Oceanography, and
  graduate student Jennifer Haas have studied the same 50-mile stretch
  of beach north, but with a completely different technique. The
  Scripps team used a mineralogical fingerprinting technique. They
  compared sand grains collected from beaches in the study area to
  grains taken from coastal bluffs, rivers, and from dredged material
  that the San Diego Regional Beach Sand Project used to replenish the
  region's disappearing beaches. 
...

http://engrish.com/

I did a quick test of trying to do increasing resolution insets so that I know where each successive SEM image is located. This works, but I did hit a limitation with Illustrator where it would not let me zoom in beyond 6400%. This is definitely a problem. It was a lot easier to place images with images successively double the resolution. Then I could just do scales by 50% repeatedly until I had the correct size. When they were not powers of 2, I had to do the scaling by eye. Not as fast or as accurate.

Not that much new in pmag-py with revision 0.4, but here it is:

• Makefile: Support for doing diffs between versions
• agfm.py: More entries from AGFM header - step and range
• agfm.py: Added --save-si, --mass (for high field chi(X)), --density

http://schwehr.org/software/pmag/

Thanks to Evelyn York for helping me through the whole process.

Yesterday was my first real time using the SIO Analytical Factilies SEM. I have sat in before while Kevin and Rodrigo ran some of their samples through, but that is not the same as going through it myself. I ran three samples through yesterday of mud that has been dried. The first task was to slice my cubes into slabs that were as flat as possible and then mount them on pedestals that have a special double stick tape.


We used the low vacuum mode on the FEI Quanta 600 E-SEM. Images were taken in both BSE and SE mode. Here is an image from the BSE mode of pyrite framboids (typical in somewhat reducing sediments).



The instrument also has an Energy Dispersive Xray spectroscopy mode (EDX) that is able to composition work. This mode gives an estimate of abundance in a 5 micron spot. You can also get quick look comparisons between two spots like this:



I did not have time to do a compositional map, but that looks like a very useful capability.

I was hoping that the tif images from the SEM would have some of the instruments status in the tag headers. Things like knowing the stage position would be really great (especially for automated visualizations), but no such luck. Here is what I get from tiffdump. tiffdump is a part of the libtiff package.

tiffdump 54_BSE_100x_5.tif 
54_BSE_100x_5.tif:
Magic: 0x4949  Version: 0x2a
Directory 0: offset 974609 (0xedf11) next 0 (0)
SubFileType (254) LONG (4) 1<0>
ImageWidth (256) SHORT (3) 1<1024>
ImageLength (257) SHORT (3) 1<948>
BitsPerSample (258) SHORT (3) 1<8>
Compression (259) SHORT (3) 1<1>
Photometric (262) SHORT (3) 1<1>
Make (271) ASCII (2) 22
Artist (315) ASCII (2) 1<\036>
StripOffsets (273) LONG (4) 1<32>
StripByteCounts (279) LONG (4) 1<970752>
SamplesPerPixel (277) SHORT (3) 1<1>
Resolution (282) RATIONAL (5) 1<200>
YResolution (283) RATIONAL (5) 1<200>
ResolutionUnit (296) SHORT (3) 1<1>
Software (305) ASCII (2) 13
33560 (0x8318) LONG (4) 1<974493>

And here is a little bit nicer display of what is in the file:

tiffinfo 54_BSE_100x_5.tif 
TIFFReadDirectory: Warning, 54_BSE_100x_5.tif: invalid TIFF directory; tags are not sorted in ascending order.
TIFFReadDirectory: Warning, 54_BSE_100x_5.tif: unknown field with tag 33560 (0x8318) encountered.
TIFF Directory at offset 0xedf11
  Subfile Type: (0 = 0x0)
  Image Width: 1024 Image Length: 948
  Resolution: 200, 200 (unitless)
  Bits/Sample: 8
  Compression Scheme: None
  Photometric Interpretation: min-is-black
  Artist: "\036"
  Make: "SIS, D-48153 Muenster"
  Samples/Pixel: 1
  Planar Configuration: single image plane
  Software: analySIS 3.2
  Tag 33560: 974493

I wanted to be able to plot both the raw AGFM and slope corrected plots on the same graph. The slope of the raw AGFM is related to the high field susceptibility and the flattened is what people are used to seeing. In a future post, I will hopefully explain how to get the normalized high field susceptibility (k(hf) or Chi_hf... need MathML!).

agfm.py is becoming a beast/garbage can. It now has --save-si in addition to --save-flattened-si. The --save-si saves the raw data in a gnuplot style ASCII format. There are some '#' comments at the top with AGFM header values and others that give the units. Then the data is space delimited columns. Hopefully that should be easy for others to parse. Here is how I generated the files to plot:

  agfm.py -d .2 --save-si --save-flattened-si bp04-1gw-s2-026c.raw 

Then I use gnuplot to see the results:

  set xlabel 'mT'
  set ylabel 'Am^2'
  plot 'bp04-1gw-s2-026c.raw.flat.si' with l, \
    'bp04-1gw-s2-026c.raw.si' with l

Which results in a plot like this: Finally, here is a little bash script that does the whole deal. It uses Lisa's hystcrunch that is in her pmag package. You should verify for yourself that the units and results are the same (or at least within fit error). Then it uses agfm.py to produce gnuplot files. Finally, I just noticed that gnuplot has svg output. Woo hoo! It is about time. I tried it and it imports into Adobe Illustrator without the problems that come with the pdf terminal. Yes, I am very excited.

#!/bin/bash
for file in *.raw; do
    echo
    echo "======== $file ======="
    agfm.py --save-flattened --save-si --save-flattened-si --agfm-format -d .2 $file 
    ls $file $file.flat
    hystcrunch -mal $file      < $file      | plotxy; mv mypost $file.ps     ; ps2pdf $file.ps
    hystcrunch -mal $file.flat < $file.flat | plotxy; mv mypost $file.flat.ps; ps2pdf $file.flat.ps
    cat << EOF | gnuplot 
set xlabel 'mT'
set ylabel 'Am^2'
set terminal gif
set output "$file.si.gif"
plot 'bp04-1gw-s2-026c.raw.flat.si' with l, 'bp04-1gw-s2-026c.raw.si' with l
set terminal pdf
set output "$file.si.pdf"
replot
set terminal svg
set output "$file.si.svg"
replot
EOF
    echo $file ... done.
done

(This is my 2nd try at this since the first draft went up in smoke
I wanted to create the Day Plot and Lisa's squareness plot for some of my hysteresis data (from the AGFM). I reverted back to gnuplot for these. From a previous post, I have created a small text file that has the hysteresis parameters generated with hystcrunch in Lisa's pmag-1.9. flattened.hyst.dat:

# flattened.hyst.dat from --dist-percent=0.2 
# filename                 sid       s        Ms        Mr*Ms   Bc     Bcr Bcrn*1e-3
bp04-1gw-s2-026.raw.flat         -.602E-37 0.599E-07 0.130E-07 0.014 0.039 0.107
bp04-1gw-s2-032.raw.flat         -.602E-37 0.106E-06 0.217E-07 0.013 0.038 0.109
bp04-1gw-s2-053.raw.flat         -.602E-37 0.225E-07 0.314E-08 0.007 0.035 0.106
...

I am trying to match Figure 13 in Tauxe et al. 2002 on micromagnetic modeling. Here is a bit of my Makefile turned into a shell script...

#!/bin/bash
# Columns in flattened.hyst.dat for awk
# Ms  = 3
# Mr  = 4
# Bc  = 5
# Bcr = 6
#
# dayplot:
egrep 'bp04-1gw' flattened.hyst.dat > hyst-1gw.dat
echo '# Bcr/Br  Mr/Ms (aka squareness)' > hyst-1gw-day.dat
awk '{print $$6/$$5,$$4/$$3}' hyst-1gw.dat >> hyst-1gw-day.dat
./dayplot.gp
#
#square-coerce:
egrep 'bp04-1gw' flattened.hyst.dat > hyst-1gw.dat
echo '# Bc(mT)  Mr/Ms (aka squareness)' > hyst-1gw-squarec.dat
awk '{print $$5*1000.,$$4/$$3}' hyst-1gw.dat >> hyst-1gw-squarec.dat
./squarec.gp

A little awk to make Bcr/Br versus Mr/Ms for the Day plot. Then the 1000 multiplier in the second awk makes Bc be in mT. I then have two gnuplot files to make the plots. Here is the dayplot.gp:

#!/usr/bin/env gnuplot
set xlabel 'Bcr/Bc'
set ylabel 'Squareness (Mr/Ms)'
set xrange [0:9.5]
set yrange [0.0:0.7]
set label 1 'SD' at 1,.55
set label 2 'PSD' at 2,0.4
set label 3 'MD' at 6,.05



set terminal pdf
set output 'bpsio2004-dayplot.pdf'
plot 'hyst-1gw-day.dat','hyst-2gw-day.dat','hyst-4gw-day.dat', \
  0.5, 0.15, 'day.l1' with l, 'day.l2' with l

It first sets the axis labels and the range of the plot. The set label puts three text labels into the plot. The coordinates are the same as the data. I don't know how to draw vertical lines with gnuplot, so I do it with the end points in file like day.l1:

1.5 0 
1.5 0.7

Then for squarec.gp:

#!/usr/bin/env gnuplot
set xlabel 'Bc (mT)'
set ylabel 'Squareness (Mr/Ms)'
set xrange [0:110]
set yrange [0.0:1.0]
set label 30 'MD' at 1,.02
set label 40 'MD+stress' at 40,.25
set label 50 'USD/SP' at 32,.35
set label 60 'vortex' at 6,.2
set label 70 'CSD' at 5,.8 # Cubic
set label 80 'flower' at 5,.70
set label 81 'flower' at 21,.4
set label 90 'complicated shapes' at 45,.65
set label 100 'uniform' at 56,.52
#
set terminal gif
set output 'bpsio2004-squarecplot.gif'
plot 'hyst-1gw-squarec.dat','hyst-2gw-squarec.dat','hyst-4gw-squarec.dat',\
  'tauxe2002-fig13.dat' with l
#
set terminal pdf
set output 'bpsio2004-squarecplot.pdf'
replot

The same basic idea, but there is a trick: tauxe2002-fig13.dat. This file contains a series of lines separated by a blank like that were digitized from Lis'a paper. I converted the pdf to a png and then used g3data (much like the old DataThief) to grab the lines.



Looks like I made it without any accidents. I hope that made some sense. matplotlib seems to be much more powerful than gnuplot, but I know gnuplot really well and am still struggling with matplotlib.

I just figured out how to quickly color a stereonet (equal area plot) that comes from Lisa's plotxy program. This is adapted from my Makefile so that it is easier to understand. The $$ and ";\" symbols required with GNU make to have embedded shell scripts gets a little confusing

#!/bin/bash
ams-eigvector-plots.py
cores=1 2 4
for core in ${CORES} ; do 
    echo $core
    plotdi < $core-ams-max.di > $core-ams-max.di.plotxy
    perl -p -e 's/symb 19 \.2/symb 0 \.075\ncolor 3\n/' $core-ams-max.di.plotxy | plotxy 
    mv mypost $core-ams-max.ps
    ps2pdf $core-ams-max.ps
done

The .py python program just creates ASCII dec inc table files that can be passed into plotdi. The big trick here is to use a perl search and replace. The plotdi output commands to pass to plotxy. Those commands are ascii. The symbol command controls the shape of the points of the stereonet. I need to replace the black dots with squares and circles. Additionally, I need to change the color to match those used for AMS eigenvector plots. Here is a little table:

  AMSEig  Eig  Color   ColorNum  Symbol    Plotxy symbol number
  Min     V3   Red      2        Circle      17 
  Int     V2   Yellow   7        Triangle     1
  Max     V1   Blue     3        Square       0

For the Maximum eigenvector, the perl line replaces the symb 19 .2 with symb 0 .075. The 0 makes the symbols be squares and the .075 makes them much smaller. Then the color 3 sets the color for the rest of the graph. This messes up the color for the error ellipse on the right hand part of the graph, but I crop that part anyway.

From the apple store...

The Power Mac G5 Quad, with a suggested retail price of $3,299 (US), includes:

• two dual-core 2.5GHz PowerPC G5 processors;
• 512MB of 533 MHz DDR2 SDRAM expandable up to 16GB;
• 250GB Serial ATA hard drive running at 7200 rpm;
• NVIDIA GeForce 6600 with 256MB of GDDR SDRAM
• three open PCI Express expansion slots: two 4-lane slots and one 8-lane slot;
• dual Gigabit Ethernet ports;
• 16x SuperDrive with double-layer support (DVD+R DL/DVD±RW/CD-RW); and ships with Mighty Mouse and Apple Keyboard.

Ways to keep it busy...

• make -j5
• 4 node rendering cluster (blender, bluemoon, maya, renderman, lightwave, povray)
• seismic data crunching
• monster volumetric rendering
• SGI, how about porting performer?
• SETI@Home
• and...

These things are going to need some serious cooling.

I remember when I got a dual 200 MHz R4K Octane MXI in 1997 for something like $50K!!!! (yikes) It really changed how I thought about computing and 3D graphics. For an order of magnitude less money, I bet this will do the same to many folks.

Hmmm... Journaling by Steven Douglas Olson on the Orielly Blogs. He talks about keeping a paper journal while you work. I keep a couple different ones at this point and go back and forth all the time.

• paper based - I have a hard bound, no lines drawing book as my paper based note book. It is easy to sketch in. I add a table of contents to the front and it is really easy to flip through. Sometimes I scan in the log book when it is really important like being used as a cruise log (Warning... that was my first attempt at a photoshop built web page and it sucks)
• blog based - This is where I try to write in a style where anyone (well, almost) can see what I am doing and is hopefully does using half-way passable English. Indexed by Google, grepped by me. This is even more valuable (to me) now that I include images when relevant.
• text file - My true garbage can of notes. I have one file per year and paste in whatever junk I want. This is not released to the world so it gets little if any actual writing, but lots of command lines and web links. Grep, grep, grep...
• readme's... - Many of my sub-projects have readme type files where very specific information goes. Sometimes this is good, other times it feels like it is diluting the power of my "one big textfile". Too many places to grep.
• revision control logs - I try to always write good logs, but it is rare that I go back to them. I need to start using some better tools! This is mostly to keep myself ready to work in groups again where reading the cvs/svn logs is a regular occurrence. Maybe I need a bridge to a news feed of log entries. That would be great! (or could really be annoying)
• wiki's - (Aurelio's fav) I think I have make 3 total entries in wiki's in my entire life. They sound good, but I think I like the other methods more, at least for now...

GMT tip from Paul Wessel:

Use pscoast -M to extract ASCII lon,lat files from the GMT cdf files.
Alternatively, get the gshhs*.b files and use the gshhs programs in
the gshhs supplement.
-pw

http://www.seacoastonline.com/news/10192005/news/68684.htm

...
  The 14-ton vessel includes a radio, a hand-held global positioning
  system, extra batteries for the GPS, two sets of nautical maps and
  personal cell phones.



  A sighting of flares was reported approximately 30 miles off
  Portsmouth Harbor Monday night, but Newlin said the Coast Guard
  searched the area with no results. There have been no reports of
  debris washed up on shore or floating on the water. 
...

I addition to wanting to know what apps are running on people's desk for how long and which is in focus, I would also like this:

http://www.oreillynet.com/pub/wlg/8137

  Screen Use Monitor: leave it running for a few days, then take a
  look at the report it generates. The app would show you with
  diagrams and animated screenshots how you use your computer's
  display; where you tend to be working most, and which bits of the
  display tend to get neglected. 

How do we actually use our computers? Once we get a handle on that, how do we improve the experience for scientists?

On a separate note... USGS to Contract Out Topo Maps

Tackling the last, great untapped wireless opportunity [gizmag.com]

  ... with all of the growth, ubiquity and connectivity, one wireless
  marketplace remains virtually untapped: the open seas. SeaMobile is
  a new venture launched in 2005 to focus on enabling wireless phones
  and computers to work at sea. SeaMobile uses scalable, flexible and
  IP/software based technology to create a user experience that is
  optimal for private yacht owners, cruise ships, ferries, offshore
  platforms, container ships and any other sea-going craft. Coverage
  can include the entire ship or select areas to enable "quiet
  zones". The technology enables people to use their own cellular
  phones and Wi-Fi enabled computers transparently while anywhere at
  sea with charges for use of the system appearing on customers'
  regular monthly bill.



  SeaMobile has established more than 200 separate roaming agreements
  with wireless carriers and systems around the world to establish
  transparent connectivity for GSM, CDMA and GPRS services. The
  SeaMobile system is "turned off" within 12 nautical miles of land to
  avoid interference with land-based wireless carrier networks.

Errr... why is "turned off" in quotes?

Just saw Woz on TV...

http://www.historychannel.com/invent/ - Win $25K.

  Has the next BIG IDEA - been your idea - but someone else got it to
  market first? Have you spent your time pondering, tinkering,
  building, or crafting the next Modern Marvel of the 21st century -
  but weren't sure what to do with it?
 
  Now you have a chance to be recognized and win a $25,000 grant for
  sharing your BIG IDEA and expressing your inventive spirit.

Nature RSS Feeds... it is not the greatest, but it is definitely a start. Where is the feed for earth and planetary sciences?

Quake aid hampered by ban on web shots [nature.com]

Open-access satellite images are revolutionizing responses to
disasters. Yet the government of Pakistan has forced aid agencies to
remove pictures of earthquake devastation from the Internet.
...
  But Anne Wright, a computer scientist at Ames Research Center who
  was involved in mapping the damage after Katrina, says the lack of
  images of Pakistan is hampering the Global Connection's efforts. "We
  haven't played a similar role with the Pakistan earthquake because
  we've had no data sources to work with. After Hurricane Katrina, the
  National Oceanic and Atmospheric Administration made their
  post-hurricane imagery publicly available. For Pakistan, the only
  available imagery we're aware of is commercial satellite
  photography."
...

http://www.kathryncramer.com/

Wow. It is actually raining reasonably hard here in la Jolla. I hope the hills turn green soon and that there is little flooding or mudsliding.

I don't really trust my agfm program correctly calcuate hysteresis parameters, so I am using Lisa's hystcrunch, but using agfm.py to flatten the curves. I have converted my Makefile into a bash script here.

 #!/bin/bash
 
 METHOD="--dist-percent=0.2"
 # flattened:
 for file in bp*.raw; do echo $file 
     agfm.py $file ${METHOD} -s ; mv $file.flat $file.flat.gnuplot
     agfm.py $file ${METHOD} -s --agfm-format 
 done
 
 # hystcrunch-flattened:
 for file in bp*.raw.flat; do 
     echo $file 
     hystcrunch -mal $file < $file | plotxy
     mv mypost $file.ps
     convert $file.ps $file.png 
 done
 
 # hystcrunch-flattened-data:
 rm -f flattened.hyst.dat
 touch flattened.hyst.dat
 echo "# flattened.hyst.dat from ${METHOD} " >> flattened.hyst.dat
 echo "# filename sid s Ms Mr*Ms Bc Bcr Bcrn*1e-3" >> flattened.hyst.dat
 for file in bp*.raw.flat; do 
     echo $file 
     echo -n "$file " >> flattened.hyst.dat 
     hystcrunch -map < $file >> flattened.hyst.dat
 done
 
 # multipage:
 rm -f multipage.{ps,pdf,plotxy}
 touch multipage.ps
 for file in bp*flat.ps; do 
     echo $file 
     cat $file >> multipage.ps 
 done
 ps2pdf multipage.ps
 open multipage.pdf

The end result of this script is a data table like this:

# flattened.hyst.dat from --dist-percent=0.2 
# filename sid s Ms Mr*Ms Bc Bcr Bcrn*1e-3
bp04-1gw-s2-026.raw.flat   -.602E-37 0.599E-07 0.130E-07 0.014 0.039 0.107
bp04-1gw-s2-032.raw.flat   -.602E-37 0.106E-06 0.217E-07 0.013 0.038 0.109
bp04-1gw-s2-053.raw.flat   -.602E-37 0.225E-07 0.314E-08 0.007 0.035 0.106
bp04-1gw-s2-062.raw.flat   -.602E-37 0.362E-07 0.395E-08 0.006 0.034 0.111
...

For a change, I used code2html with '-n' to give line nubmers. At line 3, I say that I want to use 20% of the each line segment for flattening the cure. The default of 5% was very susceptible to noise.

Then in line 10 I make a png of each hystcrunch image.

In the data section (line18) I call hystcrunch in the no plot mode (-p) so that it will pass back the parameters. I don't give anything for -l so the sid field is empty.

Finally in the section starting at line 29, I use plotxy to make a multipage pdf. plotxy has already made a postscript (ps) file for each hysteresis loop. Then I just cat these together into a large postscript file and ps2pdf magically makes a 22 page pdf. To make 8 samples per piece of paper, I print them in acrobat with 4 per page and double sided.

I am noticing that the RSS feed from GISuser.com is more like a stream of fliers for commercial packages. Not what I was hoping for. But there are still some interesting bits of information hiding in there.

http://www.gisuser.com/content/view/7160/ - Thales and LightPointe Form Alliance to Deliver Secure Outdoor Wireless Solutions

...
  This DoD mandate requires that all unclassified wireless
  transmission must be encrypted with Federal Information Processing
  Standards (FIPS) 140-2 Level 2 certified hardware utilizing the AES
  encryption algorithm. This alliance will facilitate an accelerated
  deployment of secure outdoor wireless solutions within the DoD
  agencies by enabling the use of low-cost commercial off-the-shelf
  (COTS) communication products coupled with high-bandwidth
  LightPointe optical wireless solutions for seamless, cost-effective
  connectivity. LightPointe products operate at OSI Layer 1 and are
  easily integrated with Thales' Datacryptor(R) family of encryption
  devices to meet FIPS compliance. 
...

pmag-py now has pdfposter for taking a bunch of images and laying them into a large poster. It needs lots of work, but the basics are there.

• Makefile: man pages and html versions too
• agfm.py: --{x,y}{min,max}, --ms-norm for M/Ms plots
• pdfposter.py: New program for large poster of figures

The new man pages:

• agfm
• agfmsql
• pdfposter
• units

Here is a slightly more complicated program (actually an ugly hackish one). This illustrates how to scale and write a set of png images into a large poster. It's really not that hard once you have a template. There is some trouble with worrying about image size and column wrapping that I have done a terrible job handling. Read through and try this one after looking at my previous post. When you run it, just pass in a list of png files.

I would be much smarter to watch the maximum image width for each column and jump right just enough. Also, I need to label the image file names for each one so that I don't get totally turned around. How can I denote groups/collections of input images in optparse?

#!/usr/bin/env python
__version__ = '$Revision: 1.2 $'.split()[1]
#
import math,sys
from pmag_verbosity import BOMBASTIC,VERBOSE,TRACE,TERSE,ALWAYS
from pmag_version import VERSION
import datetime
from pdflib_py import *
#
points36Inches = 2592
points60Inches = 4320
#
pointsWide = points36Inches
pointsHigh = points60Inches
#
if __name__ == "__main__":
    from optparse import OptionParser
    myparser = OptionParser(usage="%prog [options] hyst1.raw hyst2.raw hyst3.raw ...",
                            version='%prog '+__version__)
    (options,args) = myparser.parse_args()
#
    filename = 'agfm-poster.pdf'
    pdf = PDF_new()
    if -1 == PDF_open_file(pdf,filename):
        system.exit('Unable to open new pdf file '+filename)
    PDF_set_info(pdf, 'Creator', 'agfm-poster.py')
    PDF_set_info(pdf, 'Author', 'Kurt Schwehr')
    PDF_set_info(pdf, 'Title', 'AGFM of BPSIO-2004 cubes')
#
    PDF_begin_page(pdf,pointsWide,pointsHigh);
#
    font = PDF_load_font(pdf,"Helvetica-Bold","iso8859-1","")
    PDF_setfont(pdf,font, 42)
    PDF_set_text_pos(pdf,750,pointsHigh - 50)
    PDF_show(pdf,"BPSIO 2004 -- AGFM -- "+ str(datetime.date.today()));
#
    sampleHeight = 210
    sampleWidth = sampleHeight
    pointsTop = pointsHigh - 150
    y = pointsTop
    y-=36
    columnTop = y
    x = sampleHeight
#
    for file in args:
#
        print 'file = ',file
        if y < 0:
            print "Ran off the bot of the page!  Wrapping to the next column"
            y = columnTop
            x += sampleWidth
        print ' now at ',x,y
        pdfImage = PDF_load_image(pdf,"auto", file, "")
#
        if -1 != pdfImage:
            print "Putting image at ",x,y
            PDF_fit_image(pdf,pdfImage,x,y,"scale 0.2")
            PDF_close_image(pdf,pdfImage)
            y -= sampleHeight
        else: sys.exit("ERROR: unable to open " + file)
#
    PDF_end_page(pdf)
    PDF_close(pdf)
    PDF_delete(pdf)

Here is a pretty basic example for how to use pdflib to put some text in a new pdf document. Once I get the basic pdf generated, it takes some work to get everything I want in there, but it is not too hard.

#!/usr/bin/env python
import datetime
from pdflib_py import *



points32Inches = 2304
points36Inches = 2592
points60Inches = 4320
points72Inches = 5184



pointsWide = points36Inches
pointsHigh = points60Inches



if __name__ == "__main__":
    filename = 'agfm-poster.pdf'
    pdf = PDF_new()
    if -1 == PDF_open_file(pdf,filename):
        system.exit('Unable to open new pdf file '+filename)
    PDF_set_info(pdf, 'Creator', 'agfm-poster.py')
    PDF_set_info(pdf, 'Author', 'Kurt Schwehr')
    PDF_set_info(pdf, 'Title', 'AGFM of BPSIO-2004 cubes')



    PDF_begin_page(pdf,pointsWide,pointsHigh);



    font = PDF_load_font(pdf,"Helvetica-Bold","iso8859-1","")
    PDF_setfont(pdf,font, 42)
    PDF_set_text_pos(pdf,750,pointsHigh - 50)
    PDF_show(pdf,"Today's date: "+ str(datetime.date.today()));



    PDF_end_page(pdf)
    PDF_close(pdf)
    PDF_delete(pdf)

I now have agfm.py doing a lot more than before, so is time to release pmag-py 0.2. BTW, how do I keep the unittest junk out of the html pydoc file? I am only doing an import unittest!

• agfm.py: lots of cleanup, added leastsquares flattening
• agfm.py: Added lineFrom2Pts, findHystParams for Ms, Mr, Bc
• agfm.py: --show-flat, --plot-flat-file, --save-flattened
• agfm.py: --agfm-format, --dist-counts, --dist-percent
• agfm.py: Labeling of plots
• units.py: Magnetic unit conversions

http://schwehr.org/software/pmag/

I got a response back from Marine Geology saying that there is no way for them to provide an RSS feed of their title and abstracts. So to all who are reading, if you have any programming interns or such, please consider having them spend a little time writing a web scraper and then spit out an RSS feed somewhere. I'd do it myself, but I can't justify the time until I finish.

Check out Lisa's 2005 paleomagnetic book here: http://earthref.org/MAGIC/books/Tauxe/2005/

  This book is a series of 16 lectures presented as individual
  chapters. Each chapter assumes knowledge of the previous chapters so
  the chapters should be read in order. Chapter 1 begins with a review
  of the physics of magnetic fields. Maxwell's equations are
  introduced where appropriate and the magnetic units are derived from
  first principles. The conversion of units between cgs and SI
  conventions is also discussed and summarized in a handy
  table. Chapter 2 turns to the Earth's magnetic field, discussing the
  geomagnetic potential, geomagnetic elements, and the geomagnetic
  reference fields. The various magnetic poles of the Earth are also
  introduced. Chapters 3 to 7 presents the rudiments of rock magnetism
  necessary to do paleomagnetism. The most important aspect of rock
  magnetism to the working paleomagnetist is how rocks can become
  magnetized and how they can stay that way. In order to understand
  this, we begin with a discussion of what the origin of magnetism is
  in rocks, including induced and remanent magnetism, anisotropy
  energy, and mechanisms of magnetization. The magnetic properties of
  geologically important minerals are described and tabulated, as well
  as tools for their recognition. Chapter 8 provides a brief tour of
  environmental magnetism and what I have called "applied rock
  magnetism''. The nuts and bolts of how to obtain paleomagnetic
  samples and how to treat them in the laboratory is the topic of
  Chapters 9 and 10. They cover sampling strategy and routine
  laboratory procedures, including measurement and
  demagnetization. Techniques for obtaining both directional and
  paleointensity information are outlined. Chapters 11 and 12 describe
  what to do with directional data. They begin with a thorough
  discussion of Fisher [1953] statistics, describing many useful
  tests. They then illustrate how to decide whether a particular data
  set is likely to be Fisherian and Chapter 13 introduces a bootstrap
  procedure for treating data that are not. Included are many
  bootstrap tests that perform similar functions to the Fisher based
  tests that have proved so useful over the years. Magnetic tensor
  data, primarily anisotropy of magnetic susceptibility are useful in
  a number of geological applications. The acquisition and treatment
  of such data is described in Chapter 13. Traditional (Hext, 1963)
  and more modern (bootstrap) approaches are presented in some
  detail. Tests for discrimination of eigen-parameters are developed
  and illustrated with many examples. The final chapters, Chapter
  14-16, provide a whistle-stop tour of various paleomagnetic
  applications. They include a discussion of paleo-geomagnetism
  magnetostratigraphy, and apparent polar wander. 

This is mostly what I use. I sometimes also use Butler's book. PALEOMAGNETISM: Magnetic Domains to Geologic Terranes

I have been trying to deal with some hysteresis data from out Model 2900 AGFM. For some unknown reason, hystcrunch freaks out on my data, so I am wrote some code to detrend a loop. I did this by finding the slope of both the left and right tails and then removing that from the loop. Here is a first pass at it without the bells and whistles (like being able to write back out a valid file agfm file to give back to hystcrunch)

The red loop (tallest) is the original (and unflattened data) from the AGFM.. The green line is the average slope for the two tails. I used 20 points on each end, which makes 40 total for each end with both directions. The Blue '+' symbols mark the flattened curve that is created by subtracting ax+b (b=0) from the data.

Fitting lines to these groups of points is accomplished with LeastSquares module from Scientific python. It has a search algorithm that that takes a function and searches for the best fit. It is a little confusing until you have an example to copy from. My fit function is call "linear." The next trick is giving a reasonable initial condition for the search. Strange results (or long search times) can happen if you pass in really bad initial conditions on tough equations.

#!/usr/bin/env python
"""Detrending hysteresis data"""



import pylab
import agfm
from agfm import rows2cols
from Scientific.Functions.LeastSquares import leastSquaresFit



def linear(parameters,values):
    a,b = parameters
    x = values
    return (a*x + b)



distance = 20 # How many points to use on each end of the loop



sample = agfm.agfm("bp04-4gw-s1-057.raw")



data = sample.data
params=[]
initial = (1,0)  # ax+b, start with a=1 and b=0



########################################
# Right hand side start
d = data[0:distance]
param,error = leastSquaresFit(linear,initial,d)
params.append(param)



########################################
# Right hand side end
d = data[-distance:-1]
param,error = leastSquaresFit(linear,initial,d)
params.append(param)






########################################
# CENTER is the left hand side of the graph
center = len(data)/2
leftData = data[center-distance:center+distance]
param,error