/****************************************************************************/ /* */ /* GramPlot.c */ /* */ /****************************************************************************/ /****************************************************************************/ /* This program reads IAGA, GADF, Dump or WDC format data file containing */ /* magnetometer data and generates a postscript (PS) file containing a */ /* stack plot of magnetograms. */ /* */ /* Usage: */ /* GramPlot [-f] [-s] [-e | -h] [-c] [-o] [-a] |-b] [-p] [-r] [-t] [-z] */ /* [<] DataFile > EPS-File */ /* [-f Dataformat] Format of the data file. Possible values for */ /* Dataformat are IAGA (default),GADF,Dump or WDC. */ /* [-s YYMMDDHH] Time of the first record included. If missing */ /* then start of file is assumed. */ /* [-e YYMMDDHH] Time of the record not included anymore. If */ /* missing then end of file is assumed. */ /* [-h HH] Number of hours included. Either -e or -h can */ /* be specified, not both. */ /* [-c Components] List of components to be plotted (e.g. XY). */ /* Grams for different components will be plotted */ /* in separate pages. If missing then all */ /* components will be plotted. */ /* [-o 'Station list'] List of stations delimited by aposthropes. The */ /* magnetograms for different stations will be */ /* printed in the given order. In the list the */ /* stations are identified by three-letter code */ /* and separated by one space. If missing then all */ /* stations are included and stations are ordered */ /* according to geographical latitude (N to S). */ /* [-a AverageTime] Don't plot each data point but instead average */ /* values over given period (seconds). This will */ /* reduce the size of the postscript file. */ /* If missing then each data point will be plotted.*/ /* |-b BaseLineInfo] Determines the baselines used in the plots. */ /* BaseLineInfo may have the following two values: */ /* 1. H1-H2 The baseline is determined as the */ /* average field value from hour H1 to hour H2 */ /* (H2 is not included). The hour is counted */ /* from the start time of the plot. */ /* E.g -s 94120106 -b 6-7 implies that the */ /* baseline is the one hour average from */ /* 94120112 to 94120113. Note that H1 and H2 */ /* may be either positive or negative (e.g. */ /* -b -5--3 ) */ /* 2. QN The baseline is determined as the */ /* average of the quietest N hour interval in */ /* the dataplot. E.g -b Q2 implies that the */ /* program computes the max-min values for the */ /* intervals 0-2,1-3,2-4, ... and uses the */ /* interval with smallest variations to compute */ /* the baseline value. */ /* If -b is missing then the baseline is set to */ /* the average value of the field over the whole */ /* period of the plot. */ /* [-p Parameter_file] Name of the parameter file containing values of */ /* several parameters used in the plotting. If one */ /* wants to override the default values then the */ /* new values must be defined here. The structure */ /* of the parameter file may be best found by */ /* investigating the example parameter file. */ /* [-r ScaleBarLength] Length of the scale bar in nT. This defines the */ /* scale used in the plots. Physically the scale */ /* bar is about 50 mm long. If ScaleBarLength is */ /* not defined then scale is determined from */ /* FieldScale whose default value is 0.04 mm/nT. */ /* FieldScale may be defined in Parameter file. */ /* [-t TitleStr] Text used as the title of the figure. Default */ /* value is 'IMAGE Magnetometer network YYYY-MM-DD'*/ /* [-z] For FMI internal use only. This option plots the*/ /* curves using Finnish local time and also prints */ /* the texts in Finnish. The times defined in -s */ /* and -e options must be, however, in UT. */ /* DataFile Name of IAGA, GADF or WDC format data file. */ /* EPS-File Name of Encapsulated PostScript (EPS) file. */ /* */ /****************************************************************************/ /****************************************************************************/ /* Lasse Hakkinen */ /* Finnish Meteorological Institute */ /* Geophysical Research Division */ /* P.O.Box 503 */ /* FIN-00101, Helsinki, Finland */ /* e-mail: Lasse.Hakkinen@fmi.fi */ /* phone : (+358)-9-19294634 */ /* fax : (+358)-9-19294603 */ /* */ /* version 1.14 28.05.2001 */ /****************************************************************************/ /****************************************************************************/ /* Version history: */ /* */ /* 1.14 28.05.2001 Added PreviousHour and NextHour routines. These fixed */ /* a timing bug when the data did not start at full hour. */ /* 1.13 15.11.2000 Added the -z option (For FMI internal use only). */ /* 1.12 14.01.2000 Added the -t option (= Title string). */ /* 1.11 04.01.2000 Added the checking of -a option (averaging time) so */ /* that given averaging time cannot be smaller than data */ /* sampling rate. */ /* 1.10 24.11.1999 Fixed a bug in computing the coordinates of the right */ /* edge of the BoundingBox. */ /* 1.09 19.10.1999 Fixed code so that the BoundingBox is now adjusted */ /* according to the ZoomFactor. */ /* 1.08 02.07.1999 Fixed a Y2K bug in NewTime.h file which resulted in */ /* incorrect year in date strings if year >= 2000. */ /* 1.07 21.10.1998 Accommodated to the changes in Dump.h (0.1 nT to0.01 nT */ /* resolution). */ /* 1.06 12.01.1998 Added century to the default title string */ /* */ /* 1.05 04.03.1997 Added the ScaleBarLength option to command line which */ /* enables user to set the field scale. */ /* 1.04 13.11.1996 Added ability to read Dump format files. Dump format is */ /* a data format used at Nurmijarvi Observatory. */ /* 1.03 02.08.1996 Changed -b H1-H2 option so that H1 and H2 may be */ /* negative. */ /* 1.02 22.03.1996 Added the -b option for handling baseline information. */ /* 1.01 26.02.1996 Cosmetic changes in code. No apparent change in program */ /* behaviour. */ /* 1.0 30.11.1995 First official release. */ /****************************************************************************/ #include #include #include #include #include "Usage.h" #include "MagnData.h" #include "IAGA.h" #include "GADF.h" #include "Dump.h" #include "WDC.h" #include "PScript.h" #define mm_to_inch72 2.835 char *version = "1.14"; char *date = "28.05.2001"; #define HeaderLineCount 4 #define UsageLineCount 70 char *HeaderText[HeaderLineCount] = { "**************************************************************************", "This program reads in a magnetometer data file and writes out a PostScript", "file which contains magnetogram stack plots of the given data. ", "**************************************************************************", }; char *UsageText[UsageLineCount] = { " [-f] [-s] [-e | -h] [-c] [-o] [-a] [-b] [-p] [-r] [-t] [-z]", " [<] DataFile > PS_File ", " [-f Dataformat] Format of the data file. Possible values for ", " Dataformat are IAGA (default),GADF,Dump or WDC. ", " [-s YYMMDDHH] Time of the first record included. If missing ", " then start of file is assumed. ", " [-e YYMMDDHH] Time of the record not included anymore. ", " If missing then end of file is assumed. ", " [-h HH] Number of hours included. Either -e or -h can ", " be specified, not both. ", " [-c Components] List of components to be plotted (e.g. XY). ", " Grams for different components will be plotted ", " in separate pages. If missing then all ", " components will be plotted. ", " [-o 'Station list'] List of stations enclosed in quotes. The ", " magnetograms for different stations will be ", " printed in the given order. In the list the ", " stations are identified by three-letter code ", " and separated by one space. If missing then all ", " stations are included and stations are ordered ", " according to geographical latitude (N to S). ", " [-a AverageTime] Don't plot each data point but instead average ", " values over given period (seconds). This will ", " reduce the size of the postscript file. ", " If missing then each data point will be plotted.", " [-b BaseLineInfo] Determines the baselines used in the plots. ", " BaseLineInfo may have the following two values: ", " 1. H1-H2 The baseline is determined as the ", " average field value from hour H1 to hour H2 ", " (H2 is not included). The hour is counted ", " from the start time of the plot. ", " E.g -s 94120106 -b 6-7 implies that the ", " baseline is the one hour average from ", " 94120112 to 94120113. ", " H1 and H2 may be either positive or negative.", " (e.g. -b -5--3). ", " 2. QN The baseline is determined as the ", " average of the quietest N hour interval in ", " the dataplot. E.g -b Q2 implies that the ", " program computes the max-min values for the ", " intervals 0-2,1-3,2-4, ... and uses the ", " interval with smallest variations to compute ", " the baseline value. ", " If -b is missing then the baseline is set to ", " the average value of the field over the whole ", " period of the plot. ", " [-p Parameter_file] Name of the parameter file containing values of ", " several parameters used in the plotting. If one ", " wants to override the default values then the ", " new values must be defined here. The structure ", " of the parameter file may be best found by ", " investigating the example parameter file. ", " [-r ScaleBarLength] Length of the scale bar in nT. This defines the ", " scale used in the plots. Physically the scale ", " bar is about 50 mm long. If ScaleBarLength is ", " not defined then default value for field scale ", " is 0.04 mm/nT. Field scale may also be defined ", " in the parameter file. However, if -r option is ", " defined then the value given in the parameter ", " file is neglected. ", " [-t TitleString] Text to be used as the title of the figure. ", " Default title is 'IMAGE Magnetometer Network ", " YYYY-MM-DD'. If -t option is used then the value", " defined in the parameter file is neglected. ", " [-z] For FMI internal use only. Print the grams using", " Finnish local time. Also prints the texts in ", " Finnish. The times defined in -s and -e options ", " must be, however, in UT. ", " DataFile Name of IAGA, GADF, Dump or WDC data file. ", " PS_File Name of PostScript (PS) file. " }; /*==========================================================================*/ /* There are a lot of global variables here (for simplicity !). They are */ /* divided into two groups: 1. Boolean variables which control whether some */ /* parts are printed or not and 2. variables that control how different */ /* parts of the grams are printed (e.g. linewidth). */ /* Most of the variables are also given a default value which will be */ /* replaced by a value read from the parameter file (if defined). Values */ /* for some of the parameters (def value = -1) are computed during program */ /* execution. */ /*==========================================================================*/ /*------------------------*/ /* Boolean variables */ /*------------------------*/ long PrintFrame = 1; /* Print the frame enclosing the stack plot */ long PrintTitle = 1; /* Print the title string */ long PrintGridLines = 1; /* Print vertical grid lines */ long PrintTimeAxis = 1; /* Print time axis with tick marks */ long PrintLowerTicks = 1; /* Print tick marks in the lower time axis */ long PrintUpperTicks = 1; /* Print tick marks in the upper time axis */ long PrintTimeLabels = 1; /* Print time labels under the time axis */ long PrintCompChar = 1; /* Print the field component identifier */ long PrintScaleBar = 1; /* Print the field amplitude scale bar */ long PrintBaseline = 1; /* Print the baselines for each station */ long PrintBaselineValue = 1; /* Print the baseline value */ long PrintStationID = 1; /* Print the three letter Station ID. */ long PrintAverageValue = 1; /* Print the averaging value used in plots */ long PlotMode = 1; /* 1 : continuous line; 0: single dots */ long LandscapeMode = 0; /* 1 : long side horizontal 0: vertical */ /*------------------------*/ /* Other global variables */ /*------------------------*/ /* --- Title --- */ char TitleStr[100] = ""; /* Title of the magnetogram plot */ char TitleFont[20] = "Helvetica"; /* Title string font type */ float TitleFontSize = 18.0; /* Font size of title string (in points) */ float TitleOffset = 10.0; /* Distance of Title string from frame top */ /* --- Frame --- */ float FrameLeft = 40.0; /* x-coordinate of the left edge of frame */ float FrameBottom = 35.0; /* y-coordinate of the bottom edge of frame */ float FrameWidth = 140.0; /* Width of the enclosing frame */ float FrameHeight = 220.0; /* Height of the enclosing frame */ float FrameLineWidth = 0.2; /* Line width of the enclosing frame */ /* --- Data line/points --- */ float DataLineWidth = 0.05; /* Width of the dataline */ float DotRadius = 0.2; /* Radius of data points if PlotMode = 0 */ /* --- Time tickmarks --- */ long MajorTickPeriod = -1; /* Period of successive major time ticks */ long MinorTickPeriod = -1; /* Period of successive minor time ticks */ float MajorTickLength = 3.5; /* Length of the major time tick marks */ float MinorTickLength = 2.5; /* Length of the minor time tick marks */ float TickMarkWidth = 0.2; /* Width of the tick marks in time axis */ float TimeTickLeft = 0.0; /* Distance of first tick from left edge */ float TimeTickRight = 0.0; /* Distance of last tick from right edge */ /* --- Time labels --- */ char TimeLabelFont[20] = "Helvetica"; float TimeLabelFontSize = 13.0; /* Font size of time labels (in points) */ float TimeLabelOffset = 7.0; /* Distance from time labels to time axis */ float TimeTextOffset = 20.0; /* Dist of time unit text from bottom axis */ /* --- Component character --- */ char CompFont[20] = "Helvetica"; /* Field component string font type */ float CompFontSize = 16.0; /* Font size of the component string */ float CompCharOffset = 13.0; /* Distance of Comp string from Frame edge */ /* --- Field strength scale bar --- */ float FieldScale = 0.04;/* unit mm/nT */ float ScaleBarOffset = 15.0;/* Distance from scale bar to Frame edge */ float ScaleBarEndLength = 3.0;/* Length of the scale bar bottom line */ char ScaleBarFont[20] = "Helvetica"; /* Scale bar text font type */ float ScaleBarFontSize = 12.0;/* Font size of the Scale bar text */ float ScaleBarTextOffset = 3.0;/* Distance from scale bar to value text */ float ScaleBarLineWidth = 0.3;/* Line width of the scale bar */ /* --- Station identifier --- */ char StationFont[20] = "Helvetica";/* Font used in writing Station ID */ float StationFontSize = 13.0;/* Font size of the Station ID text */ float StationLabelOffset = 5.0;/* Distance of Station ID from right edge */ /* --- Baseline --- */ char BaseLineFont[20] = "Helvetica";/* Font used in writing baseline value */ float BaseLineFontSize = 9.0;/* Font size of the baseline value text */ float BaseLineLabelOffset = 2.0;/* Distance of basevalue from left edge */ float BaselineWidth = 0.05;/* Width of the horizontal baseline line */ /* --- Averaging period --- */ char AverageFont[20] = "Helvetica";/* Font used in writing averaging value */ float AverageFontSize = 9.0;/* Font size of the baseline value text */ float AverageOffset = 2.0; /* Distance of averagevalue from top edge */ /* --- Other variables --- */ float ZoomFactor = 1.0; /* Zoom entire figure by this factor */ float GridLineWidth = 0.05;/* Width of the vertical grid line */ float FirstStationOffset = 10.0;/* Extra space from top to 1st station base */ float LastStationOffset = 10.0;/* Extra space from bottom to last station */ long Finnish = 0; /* Use Finnish local time. -z option */ /*==========================================================================*/ /* Station baselines and baseline's y-coordinate may be defined in the */ /* parameter file. We must therefore define a linked list of a structure */ /* where to put these new baseline values. */ /*==========================================================================*/ struct BaselineStruct { char StationID[4]; /* Three letter station ID */ long Xbase,Ybase,Zbase; /* Baseline values for each component */ float BaseCoordinate; /* y-coordinate value of the baseline */ struct BaselineStruct *Next; /* Pointer to next new value structure */ }; typedef struct BaselineStruct *BaselineStructPtr; BaselineStructPtr BaselineList = NULL; /* List of new values read from the */ /* parameter file. */ /*==========================================================================*/ /* Function prototypes are here */ /*==========================================================================*/ long ReadParameterFile(char *ParamFileName); void FindTitleStr(char *Title,Time_sec StartTime); void GetCurrentTime(char *DateStr); long SummerTime(Time_sec T); void ComputeBoundingBox(long *left,long *bottom,long *right,long *top); void WritePSHeader(char *ComponentList); StationPtr GetStationPtr(Network_struct *NETWORK,char *StationList,long i); void SetTickPeriods(long Major, long Minor); void FindTickPeriods(Time_sec StartTime,Time_sec EndTime); void DrawFrame(void); void DrawTitle(void); void DrawTimeLabels(Time_sec StartTime, Time_sec EndTime); void DrawTimeAxis(Time_sec StartTime, Time_sec EndTime); void DrawComponentChar(char Component); void DrawScaleBar(long ScaleBarLength); float FindBaselineCoord(StationPtr Station, long Count, long Index); void GetHours(char *BaseStr,long *t0, long *t1); long GetQHour(char *BaseStr); long FindBaselineValue(StationPtr Station,char Component, Time_sec StartTime, Time_sec EndTime,char *BaseStr); void DrawBaseline(float yBase, long BaseValue); void DrawStationID(StationPtr Station, float yBase); void DrawMagnetogram(StationPtr Station,char Component,float yBase, long BaseValue,Time_sec StartTime,Time_sec EndTime,long AverageTime); void DrawAveragingValue(long AveragingValue); /*--------------------------------------------------------------------------*/ /* Skip given number of lines in the given text-file. */ /*--------------------------------------------------------------------------*/ static void SkipLines(FILE *ParamFile, long LineCount) { long i; for (i=0;i0) strcpy(ChrStr,ChrArray); /* If string length > 0 */ } /*--------------------------------------------------------------------------*/ /* Read plotting parameters from a parameter file. The structure of the */ /* parameter file is explained in the example parameter file. */ /*--------------------------------------------------------------------------*/ long ReadParameterFile(char *ParamFileName) { FILE *ParamFile; if ((ParamFile = fopen(ParamFileName,"r")) == NULL) return FileError; SkipLines(ParamFile,18); PrintFrame = GetNextLong(ParamFile); PrintTitle = GetNextLong(ParamFile); PrintGridLines = GetNextLong(ParamFile); PrintTimeAxis = GetNextLong(ParamFile); PrintLowerTicks = GetNextLong(ParamFile); PrintUpperTicks = GetNextLong(ParamFile); PrintTimeLabels = GetNextLong(ParamFile); PrintCompChar = GetNextLong(ParamFile); PrintScaleBar = GetNextLong(ParamFile); PrintBaseline = GetNextLong(ParamFile); PrintBaselineValue = GetNextLong(ParamFile); PrintStationID = GetNextLong(ParamFile); /* PrintAverageValue = GetNextLong(ParamFile);*/ PlotMode = GetNextLong(ParamFile); LandscapeMode = GetNextLong(ParamFile); SkipLines(ParamFile,7); GetNextStr(ParamFile,TitleStr); GetNextStr(ParamFile,TitleFont); TitleFontSize = GetNextFloat(ParamFile); TitleOffset = GetNextFloat(ParamFile); SkipLines(ParamFile,2); FrameLeft = GetNextFloat(ParamFile); FrameBottom = GetNextFloat(ParamFile); FrameWidth = GetNextFloat(ParamFile); FrameHeight = GetNextFloat(ParamFile); FrameLineWidth = GetNextFloat(ParamFile); SkipLines(ParamFile,2); DataLineWidth = GetNextFloat(ParamFile); DotRadius = GetNextFloat(ParamFile); SkipLines(ParamFile,2); MajorTickPeriod = GetNextLong(ParamFile); MinorTickPeriod = GetNextLong(ParamFile); MajorTickLength = GetNextFloat(ParamFile); MinorTickLength = GetNextFloat(ParamFile); TickMarkWidth = GetNextFloat(ParamFile); TimeTickLeft = GetNextFloat(ParamFile); TimeTickRight = GetNextFloat(ParamFile); SkipLines(ParamFile,2); GetNextStr(ParamFile,TimeLabelFont); TimeLabelFontSize = GetNextFloat(ParamFile); TimeLabelOffset = GetNextFloat(ParamFile); TimeTextOffset = GetNextFloat(ParamFile); SkipLines(ParamFile,2); GetNextStr(ParamFile,CompFont); CompFontSize = GetNextFloat(ParamFile); CompCharOffset = GetNextFloat(ParamFile); SkipLines(ParamFile,2); FieldScale = GetNextFloat(ParamFile); ScaleBarOffset = GetNextFloat(ParamFile); ScaleBarEndLength = GetNextFloat(ParamFile); GetNextStr(ParamFile,ScaleBarFont); ScaleBarFontSize = GetNextFloat(ParamFile); ScaleBarTextOffset = GetNextFloat(ParamFile); ScaleBarLineWidth = GetNextFloat(ParamFile); SkipLines(ParamFile,2); GetNextStr(ParamFile,StationFont); StationFontSize = GetNextFloat(ParamFile); StationLabelOffset = GetNextFloat(ParamFile); SkipLines(ParamFile,2); GetNextStr(ParamFile,BaseLineFont); BaseLineFontSize = GetNextFloat(ParamFile); BaseLineLabelOffset = GetNextFloat(ParamFile); BaselineWidth = GetNextFloat(ParamFile); /* SkipLines(ParamFile,2); GetNextStr(ParamFile,AverageFont); AverageFontSize = GetNextFloat(ParamFile); AverageOffset = GetNextFloat(ParamFile); */ SkipLines(ParamFile,2); ZoomFactor = GetNextFloat(ParamFile); GridLineWidth = GetNextFloat(ParamFile); FirstStationOffset = GetNextFloat(ParamFile); LastStationOffset = GetNextFloat(ParamFile); SkipLines(ParamFile,11); /* Add here the commands for reading the baselines */ fclose(ParamFile); return 0; } /*--------------------------------------------------------------------------*/ /* Find whether the date is before the start of IMAGE network (Sep 1,1991).*/ /* If so then it is EISCAT magnetometer cross otherwise IMAGE network. */ /*--------------------------------------------------------------------------*/ void FindTitleStr(char *Title,Time_sec StartTime) { Time_sec IMAGEstarttime; char DateStr[14]; IMAGEstarttime = StrToSecs("910901",0); if (StartTime SummerStartSec) && (T < SummerEndSec)) return 1; else return 0; } /*--------------------------------------------------------------------------*/ /* Compute the bounding box for eps-file. Bounding box is the smallest */ /* rectangle which totally encloses the figure. */ /*--------------------------------------------------------------------------*/ void ComputeBoundingBox(long *left,long *bottom,long *right,long *top) { *left = (long) (ZoomFactor*mm_to_inch72*(FrameLeft-ScaleBarOffset -ScaleBarTextOffset-0.35*ScaleBarFontSize-3.0)); *bottom = (long) (ZoomFactor*mm_to_inch72*(FrameBottom-TimeTextOffset-3.0)); *right = (long) (ZoomFactor*mm_to_inch72*(FrameLeft+FrameWidth +StationLabelOffset+StationFontSize+3.0)); *top = (long) (ZoomFactor*mm_to_inch72*(FrameBottom+FrameHeight+TitleOffset +0.35*TitleFontSize+3.0)); } /*--------------------------------------------------------------------------*/ /* Write the Title, Creation date and Bounding box into eps-file and do */ /* other PS-initializations. */ /*--------------------------------------------------------------------------*/ void WritePSHeader(char *ComponentList) { char PSTitleStr[100]; char CreationDateStr[100]; long BBleft,BBbottom,BBright,BBtop; strcpy(PSTitleStr,TitleStr); strcat(PSTitleStr,", "); strcat(PSTitleStr,ComponentList); strcat(PSTitleStr,"-component"); if (strlen(ComponentList) > 1) strcat(PSTitleStr,"s"); GetCurrentTime(CreationDateStr); ComputeBoundingBox(&BBleft,&BBbottom,&BBright,&BBtop); InitializePS(PSTitleStr,"GramPlot-program",CreationDateStr, BBleft,BBbottom,BBright,BBtop); } /*--------------------------------------------------------------------------*/ /* Find the i'th station from StationList and return a pointer to the */ /* corresponding station structure. If StationList is empty then find the */ /* i'th northernmost station in the NETWORK. */ /*--------------------------------------------------------------------------*/ StationPtr GetStationPtr(Network_struct *NETWORK,char *StationList,long i) { StationPtr Station,Station2; long Number; if (strlen(StationList) > 0) return(FindStation(NETWORK,StationList+(i-1)*4)); /* Find the i'th northernmost station in the NETWORK. */ /* This doesn't work if there are two or more stations */ /* with exactly the same latitude ! (I'll fix this later) */ Station = NETWORK->StationList; do { Number = 1; Station2 = NETWORK->StationList; while (Station2 != NULL) { if (Station->Latitude < Station2->Latitude) Number++; Station2 = Station2->Next; } if (Number != i) Station = Station->Next; } while (Number != i); return(Station); } /*--------------------------------------------------------------------------*/ /* Set the MajorTickPeriod and MinorTickPeriod. If their values are -1 */ /* then new values will be substituted, otherwise the values read from */ /* the parameter file will be used. */ /*--------------------------------------------------------------------------*/ void SetTickPeriods(long Major, long Minor) { if (MajorTickPeriod == -1) MajorTickPeriod = Major; if (MinorTickPeriod == -1) MinorTickPeriod = Minor; } /*--------------------------------------------------------------------------*/ /* Find reasonable MajorTickPeriod and MinorTickPeriod based on StartTime */ /* and EndTime. */ /*--------------------------------------------------------------------------*/ void FindTickPeriods(Time_sec StartTime,Time_sec EndTime) { long TotalHours,TotalMins; TotalMins = (EndTime-StartTime)/60; TotalHours = (EndTime-StartTime)/3600; if (TotalHours == 0) SetTickPeriods(5,1); else if (TotalHours == 1) SetTickPeriods(10,5); else if (TotalHours <= 3) SetTickPeriods(30,10); else if (TotalHours <= 6) SetTickPeriods(60,10); else if (TotalHours <= 12) SetTickPeriods(120,30); else if (TotalHours <= 24) SetTickPeriods(360,60); else if (TotalHours <= 96) SetTickPeriods(1440,240); else SetTickPeriods(1440,1440); } /*--------------------------------------------------------------------------*/ /* Draw the frame enclosing the magnetograms. */ /*--------------------------------------------------------------------------*/ void DrawFrame(void) { if (PrintFrame) { LineWidthPS(FrameLineWidth); RectPS(FrameLeft,FrameBottom,FrameWidth,FrameHeight,0.0); } } /*--------------------------------------------------------------------------*/ /* Draw the frame enclosing the magnetograms. */ /*--------------------------------------------------------------------------*/ void DrawTitle(void) { float x,y; if (PrintTitle) { x = FrameLeft+FrameWidth/2; y = FrameBottom+FrameHeight+TitleOffset; FontPS(TitleFont,TitleFontSize); TextPS(x,y,0.0,'C',TitleStr); } } /*--------------------------------------------------------------------------*/ /* Draw time labels under the time axis */ /*--------------------------------------------------------------------------*/ void DrawTimeLabels(Time_sec StartTime, Time_sec EndTime) { float x,y,dx; float TimeScale; Time_sec CurrTime; char FullDateStr[16]; char DateStr[16]; if ((! PrintTimeAxis) || (! PrintTimeLabels)) return; FontPS(TimeLabelFont,TimeLabelFontSize); TimeScale = (FrameWidth-TimeTickLeft-TimeTickRight)/(EndTime-StartTime); x = FrameLeft+TimeTickLeft; dx = 60*TimeScale*MajorTickPeriod; if (Finnish) { /* Display time in Finnish time */ StartTime += (2 + SummerTime(StartTime))*3600; EndTime += (2 + SummerTime(StartTime))*3600; } CurrTime = StartTime; do { SecsToStr(CurrTime,FullDateStr); strcpy(DateStr,""); if ((EndTime-StartTime) <= 86400L) {/* Print hours, not day numbers */ strncat(DateStr,FullDateStr+6,2); if (MajorTickPeriod < 60) { /* Print minutes */ strcat(DateStr,":"); strncat(DateStr,FullDateStr+8,2); } TextPS(x,FrameBottom-TimeLabelOffset,0,'C',DateStr); } else { /* Print day numbers, not hours */ strncat(DateStr,FullDateStr+4,2); /* If MajorTickPeriod is one day print day */ /* number between the tick marks */ if (MajorTickPeriod == 1440) { if (CurrTime < EndTime) TextPS(x+dx/2,FrameBottom-TimeLabelOffset,0,'C',DateStr); } else TextPS(x,FrameBottom-TimeLabelOffset,0,'C',DateStr); } x += dx; CurrTime += 60*MajorTickPeriod; } while (CurrTime <= EndTime); /* Print the text indicating time unit */ x = FrameLeft+FrameWidth/2; y = FrameBottom-TimeTextOffset; FontPS(TimeLabelFont,TimeLabelFontSize); if ((EndTime-StartTime) <= 86400L) { /* Print hour text */ if (Finnish) TextPS(x,y,0.0,'C',"Kellonaika Suomen aikaa"); else TextPS(x,y,0.0,'C',"Hour (UT)"); } else { /* Print day text */ if (Finnish) TextPS(x,y,0.0,'C',"Kuukauden pŠivŠ"); else TextPS(x,y,0.0,'C',"Day of month"); } } /*--------------------------------------------------------------------------*/ /* Draw the time axis, tick marks and vertical grid lines */ /*--------------------------------------------------------------------------*/ void DrawTimeAxis(Time_sec StartTime, Time_sec EndTime) { float x,dx; float TimeScale; Time_sec CurrTime; if (!PrintTimeAxis) return; TimeScale = (FrameWidth-TimeTickLeft-TimeTickRight)/(EndTime-StartTime); CurrTime = StartTime; x = FrameLeft+TimeTickLeft; dx = 60*TimeScale*MinorTickPeriod; do { if ((CurrTime-StartTime) % (60*MajorTickPeriod) == 0) { /* Major */ if (PrintGridLines) { LineWidthPS(GridLineWidth); LinePS(x,FrameBottom,x,FrameBottom+FrameHeight); } LineWidthPS(TickMarkWidth); if (PrintLowerTicks) LinePS(x,FrameBottom,x,FrameBottom+MajorTickLength); if (PrintUpperTicks) LinePS(x,FrameBottom+FrameHeight,x, FrameBottom+FrameHeight-MajorTickLength); } else { /* Minor */ if (PrintLowerTicks) LinePS(x,FrameBottom,x,FrameBottom+MinorTickLength); if (PrintUpperTicks) LinePS(x,FrameBottom+FrameHeight,x, FrameBottom+FrameHeight-MinorTickLength); } x += dx; CurrTime += 60*MinorTickPeriod; } while (CurrTime <= EndTime); } /*--------------------------------------------------------------------------*/ /* Draw text displaying the current component in the left part of the plot */ /*--------------------------------------------------------------------------*/ void DrawComponentChar(char Component) { float x,y; char CompStr[20]; if (!PrintCompChar) return; if (Finnish) strcpy(CompStr,"?-komponentti"); else strcpy(CompStr,"?-COMPONENT"); FontPS(CompFont,CompFontSize); x = FrameLeft-CompCharOffset; y = FrameBottom+0.4*FrameHeight; CompStr[0] = Component; TextPS(x,y,90.0,'R',CompStr); } /*--------------------------------------------------------------------------*/ /* Draw text displaying the current averaging value used in the data plots */ /*--------------------------------------------------------------------------*/ void DrawAveragingValue(long AveragingValue) { float x,y; char AverageStr[100]; if (!PrintAverageValue) return; FontPS(AverageFont,AverageFontSize); x = FrameLeft+FrameWidth; y = FrameBottom+FrameHeight+AverageOffset; if (AveragingValue < 60) { if (Finnish) sprintf(AverageStr,"%d sekunnin keskiarvot",AveragingValue); else sprintf(AverageStr,"%d second averages",AveragingValue); } else { if (Finnish) sprintf(AverageStr,"%d minuutin keskiarvot",AveragingValue/60); else sprintf(AverageStr,"%d minute averages",AveragingValue/60); } TextPS(x,y,0.0,'R',AverageStr); } /*--------------------------------------------------------------------------*/ /* Draw the scale bar displaying the magnitude of the field in the left */ /* part of the plot */ /*--------------------------------------------------------------------------*/ void DrawScaleBar(long ScaleBarLength) { float x,y,dx,dy; char BarLengthStr[12]; if (!PrintScaleBar) return; /* Determine the length of the scale bar based on FieldScale */ if (ScaleBarLength == 0) { if (FieldScale <= 0.05) ScaleBarLength = 1000; else if (FieldScale <= 0.1) ScaleBarLength = 500; else if (FieldScale <= 0.2) ScaleBarLength = 250; else if (FieldScale <= 0.5) ScaleBarLength = 100; else if (FieldScale <= 1.0) ScaleBarLength = 50; else ScaleBarLength = 10; } /* The scale bar vector */ LineWidthPS(ScaleBarLineWidth); dy = ScaleBarLength*FieldScale; x = FrameLeft-ScaleBarOffset; y = FrameBottom+0.6*FrameHeight; LinePS(x,y,x,y+dy); /* End bars */ dx = ScaleBarEndLength; LinePS(x-dx/2,y,x+dx/2,y); y = y+dy; LinePS(x-dx/2,y,x+dx/2,y); /* The scale bar length text */ x = x-ScaleBarTextOffset; y = y-dy/2; sprintf(BarLengthStr,"%d nT",ScaleBarLength); FontPS(ScaleBarFont,ScaleBarFontSize); TextPS(x,y,90.0,'C',BarLengthStr); } /*--------------------------------------------------------------------------*/ /* Find the y-coordinate of the baseline for given station. If the baseline*/ /* is not defined in the parameter file then compute a default value. */ /*--------------------------------------------------------------------------*/ float FindBaselineCoord(StationPtr Station, long Count, long Index) { float yBase; float StationDist; /* Distance between two successive stations */ BaselineStructPtr p; /* Compute the default y-value for the baseline */ StationDist = (FrameHeight-FirstStationOffset-LastStationOffset)/Count; yBase = FrameBottom+FrameHeight-FirstStationOffset -StationDist/2-(Index-1)*StationDist; /* Check, if the baseline value is defined in the BaselineList */ p = BaselineList; while (p != NULL) { if (strncmp(p->StationID,Station->StationID,3) != 0) { if (p->BaseCoordinate > 0) yBase = p->BaseCoordinate; } p = p->Next; } return (yBase); } /*--------------------------------------------------------------------------*/ /* Get the start and end hours from the BaseStr. The format of BaseStr is */ /* HH-HH. */ /*--------------------------------------------------------------------------*/ void GetHours(char *BaseStr,long *t0, long *t1) { long i = 0; long t = 0; long negative = 0; if (BaseStr[0] == '-') { negative = 1; i++; } while ((BaseStr[i] != '-') && (i < strlen(BaseStr))) { t = 10*t + (BaseStr[i]-48); i++; } if (negative) t = -t; *t0 = t; t = 0; negative = 0; if (BaseStr[++i] == '-') { negative = 1; i++; } while (i < strlen(BaseStr)) t = 10*t + (BaseStr[i++]-48); if (negative) t = -t; *t1 = t; } /*--------------------------------------------------------------------------*/ /* Get the hour count from the BaseStr. The format of the string is QHH. */ /*--------------------------------------------------------------------------*/ long GetQHour(char *BaseStr) { long i = 0; long t = 0; while (++i < strlen(BaseStr)) t = 10*t + (BaseStr[i]-48); return t; } /*--------------------------------------------------------------------------*/ /* Find the baseline value for given station and given field component. */ /* If the baseline value is not defined in the parameter file then compute */ /* the value as the average value over the whole dataplot period. */ /*--------------------------------------------------------------------------*/ long FindBaselineValue(StationPtr Station,char Component, Time_sec StartTime, Time_sec EndTime, char *BaseStr) { BaselineStructPtr p; long Base = -1; long t,t0,t1; long Hours,Max,Min; long MaxMin = 1000000; long Good; /* Check, if the baseline value is defined in the BaselineList */ p = BaselineList; while (p != NULL) { if (strncmp(p->StationID,Station->StationID,3) != 0) { switch (Component) { case 'H' : case 'X' : if (p->Xbase != 999999) Base = p->Xbase; break; case 'D' : case 'Y' : if (p->Ybase != 999999) Base = p->Ybase; break; case 'Z' : if (p->Zbase != 999999) Base = p->Zbase; break; } } p = p->Next; } /* If the baseline was not defined in the BaselineList then */ /* check the BaseStr (-b option) and compute the baseline */ /* accordingly. */ if (strlen(BaseStr) != 0) { if (BaseStr[0] == 'Q') { /* Quiet period average */ Hours = GetQHour(BaseStr); t0 = StartTime; for (t = StartTime;t Max-Min) && (Good > 80)) { MaxMin = Max-Min; t0 = t; } } Base = ComputeAverage(Station,Component,t0,3600*Hours); } else { /* Average over specified hours */ GetHours(BaseStr,&t0,&t1); Base = ComputeAverage(Station,Component,StartTime+3600*t0, 3600*(t1-t0)); if (Base == MissingValue) Base = -1; } } /* If the baseline was not defined in the BaselineList nor in */ /* the BaseStr then compute the baseline as average value over */ /* the whole period of the dataplot. */ if (Base == -1) { Base = ComputeAverage(Station,Component,StartTime,EndTime-StartTime); if (Base != MissingValue) Base = 10*((Base+5)/10); /* Round to 1 nT */ } return (Base); } /*--------------------------------------------------------------------------*/ /* Draw a horizontal dashed line indicating the field baseline and the */ /* baseline value. */ /*--------------------------------------------------------------------------*/ void DrawBaseline(float yBase, long BaseValue) { float x,y; char BaseStr[12]; if (!PrintBaseline) return; LineWidthPS(BaselineWidth); LineTypePS(1,8); LinePS(FrameLeft,yBase,FrameLeft+FrameWidth,yBase); if (!PrintBaselineValue) return; FontPS(BaseLineFont,BaseLineFontSize); x = FrameLeft-BaseLineLabelOffset; y = yBase; sprintf(BaseStr,"%d",BaseValue/10); TextPS(x,y,90.0,'C',BaseStr); } /*--------------------------------------------------------------------------*/ /* Draw the three letter Station code. */ /*--------------------------------------------------------------------------*/ void DrawStationID(StationPtr Station, float yBase) { float x,y; if (!PrintStationID) return; FontPS(StationFont,StationFontSize); x = FrameLeft+FrameWidth+StationLabelOffset; y = yBase-0.15*StationFontSize; TextPS(x,y,0.0,'L',Station->StationID); } /*--------------------------------------------------------------------------*/ /* Here is the most important routine, the drawing of a single magnetogram.*/ /* The y-coordinate of the baseline value is at yBase. If AverageTime is */ /* larger than the sampling rate then compute averages and plot them. */ /*--------------------------------------------------------------------------*/ void DrawMagnetogram(StationPtr Station,char Component,float yBase, long BaseValue,Time_sec StartTime,Time_sec EndTime,long AverageTime) { float x,y,dx; float TimeScale; Time_sec CurrTime = StartTime; long first = 1; long DotCount = 0; long *DataPtr; long Value; LineWidthPS(DataLineWidth); LineTypePS(0,0); TimeScale = (FrameWidth-TimeTickLeft-TimeTickRight)/(EndTime-StartTime); x = FrameLeft+TimeTickLeft; dx = AverageTime*TimeScale; if (StartTime < Station->StartTime) { x += (Station->StartTime-StartTime)*TimeScale; CurrTime = Station->StartTime; } if (EndTime > Station->EndTime) EndTime = Station->EndTime; /* Find the pointer to the data point corresponding to CurrTime */ switch (Component) { case 'H' : case 'X' : DataPtr = Station->X; break; case 'D' : case 'Y' : DataPtr = Station->Y; break; case 'Z' : DataPtr = Station->Z; break; } DataPtr += (CurrTime-Station->StartTime)/(Station->TimeStep); while (CurrTime < EndTime) { if (Station->TimeStep == AverageTime) Value = *DataPtr++; else Value = ComputeAverage(Station,Component,CurrTime,AverageTime); if (PlotMode == 1) { /* Draw continuous line between data points */ if (Value != MissingValue) { y = yBase+FieldScale*(Value-BaseValue)/10; if (first) { MoveToPS(x,y); first = 0; } else { LineToPS(x,y); DotCount++; if (DotCount == 300) { /* This prevents stack overflow */ ClosePathPS(); /* in many PostScript printers */ MoveToPS(x,y); DotCount = 0; } } } else { if (!first) { ClosePathPS(); first = 1; } } } else { /* Draw every data point as single dot */ if (Value != MissingValue) { y = yBase+FieldScale*(Value-BaseValue)/10; FillCirclePS(x,y,DotRadius,0.0,0); } } x += dx; CurrTime += AverageTime; } if (!first) ClosePathPS(); } /*--------------------------------------------------------------------------*/ /* Round the given time to previous full hour or next full hour. If the */ /* time is already at full hour then don't do anything. */ /*--------------------------------------------------------------------------*/ Time_sec PreviousHour(Time_sec T) { Time_struct MyTime; SecsToTm(T,&MyTime); MyTime.tm_min = 0; MyTime.tm_sec = 0; return (TmToSecs(&MyTime)); } Time_sec NextHour(Time_sec T) { Time_struct MyTime; SecsToTm(T,&MyTime); if ((MyTime.tm_min == 0) && (MyTime.tm_sec == 0)) return (T); else return (PreviousHour(T+3600)); } /*--------------------------------------------------------------------------*/ /* The main procedure */ /*--------------------------------------------------------------------------*/ int main(int argc, char *argv[]) { long params; /* Dummy index variable */ long status = 0; /* Result code for file reading */ long FileCount = 0; /* Number of files. Must be 0 or 1 */ Network_struct NETWORK; /* Read all data into this structure */ Time_sec StartTime,EndTime; /* Start and End time in seconds */ long CompIndex,StatIndex; /* Dummy indices */ long CompFound; /* Flag for checking magnetic components*/ char Component; /* Current field component */ StationPtr Station; /* Current station */ float yBase; /* Y-coordinate of the baseline line */ long BaseValue; /* Value of the field baseline (0.1 nT) */ /* Parameters read from the command line */ long HourCount = 0; long AverageTime = 0; long ScaleBarLength = 0; /* Length of the scale bar in nT's */ char FormatStr[10] = "IAGA"; char FileName[100] = ""; char StartTimeStr[14] = ""; char EndTimeStr[14] = ""; char StationList[400] = ""; char ParamFileName[100] = ""; char ComponentList[5] = ""; char BaselineStr[10] = ""; char HeaderStr[100] = ""; /*==========================*/ /* Analyse the command line */ /*==========================*/ if (argc == 1) { /* No arguments, show the usage text */ PrintUsage(argv[0],0,HeaderLineCount,UsageLineCount); return OK; } for (params = 1; params < argc; params++) { if (*argv[params] != '-') { strcpy(FileName,argv[params]); FileCount++; } else { switch (*(argv[params]+1)) { case 'f' : strcpy(FormatStr,argv[++params]); break; case 's' : strcpy(StartTimeStr,argv[++params]); break; case 'e' : strcpy(EndTimeStr,argv[++params]); break; case 'h' : HourCount = atol(argv[++params]); break; case 'a' : AverageTime = atol(argv[++params]); break; case 'r' : ScaleBarLength = atol(argv[++params]); break; case 'c' : strcpy(ComponentList,argv[++params]); break; case 'o' : strcpy(StationList,argv[++params]); break; case 'b' : strcpy(BaselineStr,argv[++params]); break; case 'p' : strcpy(ParamFileName,argv[++params]); break; case 't' : strcpy(HeaderStr,argv[++params]); break; case 'z' : Finnish = 1; break; default : fprintf(stderr,"\n### %s: \"%s\" is not an option.\n\n", argv[0], argv[params]); PrintUsage(argv[0],1,HeaderLineCount,UsageLineCount); return FAIL; break; } } } /*====================================*/ /* Check the values of the parameters */ /*====================================*/ if (FileCount > 1) { PrintUsage(argv[0],1,HeaderLineCount,UsageLineCount); return FAIL; } if ((HourCount != 0) && (strlen(EndTimeStr) > 0)) { PrintUsage(argv[0],1,HeaderLineCount,UsageLineCount); return FAIL; } if ((HourCount != 0) && (strlen(StartTimeStr) == 0)) { PrintUsage(argv[0],1,HeaderLineCount,UsageLineCount); return FAIL; } if (HourCount != 0) { SecsToStr(StrToSecs(StartTimeStr,0)+3600*HourCount,EndTimeStr); } /*==========================================*/ /* Read the data from data file into memory */ /*==========================================*/ if (!strcmp(FormatStr,"IAGA")) { status = ReadIAGA(FileName,&NETWORK,NULL,NULL,StationList); } else if (!strcmp(FormatStr,"GADF")) { status = ReadGADF(FileName,&NETWORK,NULL,NULL,StationList); } else if (!strcmp(FormatStr,"Dump")) { status = ReadDump(FileName,&NETWORK,NULL,NULL,StationList); } else if (!strcmp(FormatStr,"WDC")) { status = ReadWDC(FileName,&NETWORK,NULL,NULL,StationList); } else fprintf(stderr,"Illegal format type: %s\n",FormatStr); if (status != 0) { if (status == FileError) fprintf(stderr,"Error in reading data file\n"); if (status == OutOfMemory) fprintf(stderr,"Out of memory while reading data file\n"); if (status == FileFormatError) fprintf(stderr,"Wrong file format in input file\n"); return FAIL; } if (!strcmp(FormatStr,"Dump")) /* change units from 0.01 nT */ ChangeUnitsDOWN(&NETWORK); /* to 0.1 nT */ /*==========================================================*/ /* Determine the start,end and averaging times for the plot */ /*==========================================================*/ if (strlen(StartTimeStr) == 0) { StartTime = NETWORK.StationList->StartTime; StartTime = PreviousHour(StartTime); /* Round to the previous full hour */ } else StartTime = StrToSecs(StartTimeStr,0); if (strlen(EndTimeStr) == 0) { EndTime = NETWORK.StationList->EndTime; EndTime = NextHour(EndTime); /* Round to next full hour */ } else EndTime = StrToSecs(EndTimeStr,0); if (AverageTime == 0) AverageTime = NETWORK.StationList->TimeStep; if (AverageTime < NETWORK.StationList->TimeStep) { fprintf(stderr,"Given average time (-a) %d s is smaller that sampling rate %d s \n", AverageTime,NETWORK.StationList->TimeStep); fprintf(stderr,"No data is written into postscript file \n"); return FAIL; } /*================================================*/ /* Check that data exist for the given components */ /*================================================*/ if (strlen(ComponentList) == 0) { strcpy(ComponentList,NETWORK.StationList->Components); } else { for (CompIndex = 0;CompIndex < strlen(ComponentList);CompIndex++) { CompFound = 0; Component = ComponentList[CompIndex]; if (Component == NETWORK.StationList->Components[0]) CompFound = 1; if (Component == NETWORK.StationList->Components[1]) CompFound = 1; if (Component == NETWORK.StationList->Components[2]) CompFound = 1; if (CompFound == 0) { fprintf(stderr,"No data for component %c\n",Component); return FAIL; } } } FindTitleStr(TitleStr,StartTime); /*==============================================*/ /* Read the parameter file if necessary */ /*==============================================*/ if (strlen(ParamFileName) != 0) { status = ReadParameterFile(ParamFileName); if (status == FileError) { fprintf(stderr,"Error in reading parameter file :%s\n", ParamFileName); return FAIL; } } /*======================================================*/ /* If HeaderStr is defined by -t option in the command */ /* line the use it as the figure title instead of what */ /* ever is defined in the parameter file. */ /*======================================================*/ if (strlen(HeaderStr) > 0) strcpy(TitleStr,HeaderStr); /*===============================================*/ /* If ScaleBarLength is defined in the command */ /* line then determine the FieldScale. */ /*===============================================*/ if (ScaleBarLength != 0) FieldScale = 50.0/ScaleBarLength; /* mm/nT */ /*===============================================*/ /* Find the Major and Minor TickPeriods based on */ /* StartTime and EndTime. */ /*===============================================*/ FindTickPeriods(StartTime,EndTime); /*==============================================*/ /* Write the Title, Creator and Creation date */ /* and initialize the postscript file */ /*==============================================*/ WritePSHeader(ComponentList); /*===============================================*/ /* Plot the stations in the order of increasing */ /* geomagnetic latitude instead of geographic */ /* latitude. */ /*===============================================*/ /* if (GeoMagn != 0) ReOrderStations(&NETWORK); */ /*============================================*/ /* Go through all components and all stations */ /*============================================*/ for (CompIndex = 0;CompIndex < strlen(ComponentList);CompIndex++) { Component = ComponentList[CompIndex]; NewPagePS(CompIndex+1); if (LandscapeMode) LandscapePS(); ZoomPS(ZoomFactor); DrawFrame(); DrawTitle(); DrawTimeAxis(StartTime,EndTime); DrawTimeLabels(StartTime,EndTime); DrawComponentChar(Component); DrawScaleBar(ScaleBarLength); DrawAveragingValue(AverageTime); for (StatIndex=1;StatIndex <= NETWORK.StationCount;StatIndex++) { Station = GetStationPtr(&NETWORK,StationList,StatIndex); yBase = FindBaselineCoord(Station,NETWORK.StationCount,StatIndex); BaseValue= FindBaselineValue(Station,Component,StartTime,EndTime, BaselineStr); DrawBaseline(yBase,BaseValue); DrawStationID(Station,yBase); DrawMagnetogram(Station,Component,yBase,BaseValue, StartTime,EndTime,AverageTime); } ClosePagePS(); } ClosePS(); FreeNetwork(&NETWORK); return OK; }