/****************************************************************************/ /* */ /* StationPlot.c */ /* */ /****************************************************************************/ /****************************************************************************/ /* This program reads an IAGA, GADF, Dump, WDC or Text format data file */ /* containing magnetometer data and generates a postscript (PS) file which */ /* displays magnetograms of all components in the same page. Grams for each */ /* station are plotted in separate pages. */ /* */ /* Usage: */ /* StationPlot [-f] [-s] [-e | -h] [-o] [-a] [-b] [-c] [-r] [-p] [-z] [-t] */ /* [<] DataFile > PS-File */ /* [-f Dataformat] Format of the data file. Possible values for */ /* Dataformat are IAGA (default),GADF,Dump, WDC or */ /* TEXT. */ /* [-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. */ /* [-o 'Station list'] List of stations delimited by aposthropes. In */ /* the list the stations are identified by */ /* three-letter code and separated by one space. */ /* If missing then all stations are included. */ /* [-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. */ /* [-c Components] Plot the selected components. Possible values */ /* are X,Y,Z,H,D,F. Default value is 'XYZ'. */ /* [-r FieldAmpl] Length of vertical axis in nT. This defines the */ /* are X,Y,Z,H,D,F. Default value is 'XYZ'. */ /* [-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. */ /* [-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. */ /* [-t TitleStr] Text used as the title of the figure. Default */ /* value is 'StationName YYYY-MM-DD' */ /* DataFile Name of IAGA, GADF, Dump or WDC data file. */ /* PS-File Name of PostScript (PS) file. */ /* */ /****************************************************************************/ /****************************************************************************/ /* Lasse Hakkinen */ /* Finnish Meteorological Institute */ /* Department of Geophysics */ /* P.O.Box 503 */ /* FIN-00101, Helsinki, Finland */ /* e-mail: Lasse.Hakkinen@fmi.fi */ /* phone : (+358)-9-19294634 */ /* fax : (+358)-9-19294603 */ /* */ /* version 1.08 24.11.2004 */ /****************************************************************************/ /****************************************************************************/ /* Version history: */ /* */ /* 1.08 01.04.2004 Added HDF plotting + modified -c option to allow */ /* arbitrary number of components to be plotted. */ /* 1.07 08.12.2003 Added -z option (Finnish time + texts). */ /* 1.06 18.11.2003 Added possibility to read Text-format files. */ /* 1.05 28.05.2001 Added PreviousHour and NextHour routines. These fixed */ /* a timing bug when the data did not start at full hour. */ /* 1.04 29.05.2000 Fixed a bug due to which the -r option did not work. */ /* 1.03 19.10.1999 Fixed code so that the BoundingBox is now adjusted */ /* according to the ZoomFactor. */ /* 1.02 09.09.1999 Fixed a Y2K bug in NewTime.h file which resulted in */ /* incorrect year in date strings if year >= 2000. */ /* 1.01 23.10.1998 Adjusted for the change of field resolution (from 0.1 nT*/ /* to 0.01 nT) in the Dump.h file. */ /* 1.0 12.01.1998 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 "TextFormat.h" #include "PScript.h" #define mm_to_inch72 2.835 char *version = "1.08"; char *date = "24.11.2005"; #define HeaderLineCount 5 #define UsageLineCount 63 char *HeaderText[HeaderLineCount] = { "**************************************************************************", " This program produces magnetogram plots of given magnetic components of ", " given stations and for specified time interval (one station per page). ", " The program writes the grams in postscript format into stdout. ", "**************************************************************************", }; char *UsageText[UsageLineCount] = { " [-f] [-s] [-e | -h] [-o] [-a] [-b] [-c] [-r] [-p] [-z] ", " [<] DataFile > PS_File ", " [-f Dataformat] Format of the data file. Possible values for ", " Dataformat are IAGA (default),GADF,Dump, WDC or ", " TEXT. ", " [-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. ", " [-o 'Station list'] List of stations enclosed in quotes. ", " In the list the stations are identified by ", " three-letter code and separated by one space. ", " If missing then all stations are included. ", " [-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. ", " [-c Components] Components to be plotted. Possible values are ", " X,Y,Z,H,D,F. Default value is 'XYZ'. ", " [-r FieldAmplitude] Length of vertical axis in nT. This defines the ", " scale used in the plots. If FieldAmplitude is ", " not defined then the grams will be autoscaled ", " so that all components will totally fit inside ", " their boxes. ", " [-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. ", " [-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. ", " [-t TitleString] Text to be used as the title of the figure. ", " Default title is 'StationCode YYYY-MM-DD H1-H2. ", " If -t option is used then the value defined in ", " the parameter file is neglected. ", " 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 gram */ long PrintTitle = 1; /* Print the title string */ long PrintVerGridLines = 1; /* Print vertical grid lines */ long PrintHorGridLines = 1; /* Print horizontal 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 PrintTimeLabelsXY = 1; /* Print time labels in upper plots */ long PrintTimeLabelsZ = 1; /* Print time labels in the bottom plot */ long PrintCompChar = 1; /* Print the field component identifier */ long PrintBaseline = 0; /* Print the baselines for each component */ long PrintFieldValues = 1; /* Print the field values in the ver. axis */ 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 = 8.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-coord. of the bottom edge of Z-frame */ float FrameWidth = 140.0; /* Width of the magnetogram plot frame */ float FrameHeight = 70.0; /* Height of the frame of single component */ float FrameDistance = 8.0; /* Distance of frames of diff. components */ float FrameLineWidth = 0.2; /* Line width of the enclosing frame */ /* --- Data line/points --- */ float DataLineWidth = 0.1; /* 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.0; /* 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 = 12.0; /* Font size of time labels (in points) */ float TimeLabelOffset = 5.0; /* Distance of time labels from time axis */ float TimeTextOffset = 12.0; /* Dist of time unit text from bottom axis */ /* --- Field axis tickmarks --- */ long MajorFieldPeriod = -1; /* Period of successive major field ticks */ long MinorFieldPeriod = -1; /* Period of successive minor field ticks */ float MajorFieldLength = 3.5; /* Length of the major field axis tick mark */ float MinorFieldLength = 2.0; /* Length of the minor field axis tick mark */ float FieldTickMarkWidth = 0.2;/* Width of the tick marks in field axis */ float FieldTickTop = 0.0; /* Distance of topmost tick from top edge */ float FieldTickBottom = 0.0; /* Distance of lowest tick from bottom edge */ /* --- Field axis labels --- */ char FieldLabelFont[20] = "Helvetica"; float FieldLabelFontSize = 10.0;/* Font size of field value labels */ float FieldLabelOffset = 3.0;/* Distance of field values from field axis */ /* --- Component character --- */ char CompFont[20] = "Helvetica"; /* Field component string font type */ float CompFontSize = 16.0; /* Font size of the component string */ float CompCharOffset = 6.0; /* Distance of Comp string from Frame edge */ /* --- Baseline --- */ 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 and hor. grid line */ long GridLineType = 5; /* Line type of grid lines. see Pscript.h */ float GridLinePeriod = 3.0; /* Length of one period in dashed grid line */ long ComponentCount; /* Number of components (= graphs) per page */ long Finnish = 0; /* Use Finnish local time. -z option */ /* --- FMI logo variables --- */ char FMIlogoFont[20] = "Helvetica"; float FMIlogoFontSize = 10.0; float FMIlogoSize = 10.0; float FMIlogoOffset = 15.0; float FMIlogoTextOffset = 2.0; /*==========================================================================*/ /* Station baselines 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 */ 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 */ /*==========================================================================*/ static void SkipLines(FILE *ParamFile, long LineCount); static long GetNextLong(FILE *ParamFile); static float GetNextFloat(FILE *ParamFile); static void GetNextStr(FILE *ParamFile,char *ChrStr); static long ReadLine(FILE *ParamFile,char *Line); long ReadParameterFile(char *ParamFileName); long SummerTime(Time_sec T); void WritePSHeader(char *StationName, char *ComponentStr); void SetTickPeriods(long Major, long Minor); void FindTickPeriods(Time_sec StartTime,Time_sec EndTime); void SetFieldTickPeriods(long Major, long Minor); void DrawTitle(void); void DrawFrame(float yBottom); void FindTitleStr(StationPtr Station, char *Title, Time_sec StartTime,Time_sec EndTime); void GetCurrentTime(char *DateStr); void ComputeBoundingBox(long *left,long *bottom,long *right,long *top, long FrameCount); void DrawTimeLabels(float yBottom, Time_sec StartTime, Time_sec EndTime); void DrawTimeUnits(Time_sec StartTime, Time_sec EndTime); void DrawTimeAxis(float yBottom, Time_sec StartTime, Time_sec EndTime); void DrawFieldAxis(float yBottom, long BaseValue, long Amplitude); void DrawDeclinationAxis(float yBottom, long BaseValue, long Amplitude); void DrawComponentChar(float yBottom, char Component); void DrawAveragingValue(long AveragingValue); void GetHours(char *BaseStr,long *t0, long *t1); long GetQHour(char *BaseStr); void DrawBaseline(float yBase); long RoundFieldAmplitude(long Amplitude); long RoundBaseValue(long BaseValue,long Amplitude); void DrawMagnetogram(StationPtr Station,char Component,float yBase, long BaseValue,long Amplitude,Time_sec StartTime, Time_sec EndTime,long AverageTime); long FindBaselineValue(StationPtr Station,char Component, Time_sec StartTime,Time_sec EndTime, char *BaseStr); long FindBaseline(StationPtr Station,Time_sec StartTime,Time_sec EndTime,char Component,long MaxAmpl); /*--------------------------------------------------------------------------*/ /* 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 one line from the parameter file */ /*--------------------------------------------------------------------------*/ static long ReadLine(FILE *ParamFile,char *Line) { int c; while (((c = getc(ParamFile)) != EOF) && (c != '\n')) *Line++ = c; *Line = '\0'; if (c == '\n') return 1; else return 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) { BaselineStructPtr p,q; char Line[100]; FILE *ParamFile; if ((ParamFile = fopen(ParamFileName,"r")) == NULL) return FileError; SkipLines(ParamFile,18); PrintFrame = GetNextLong(ParamFile); PrintTitle = GetNextLong(ParamFile); PrintVerGridLines = GetNextLong(ParamFile); PrintHorGridLines = GetNextLong(ParamFile); PrintTimeAxis = GetNextLong(ParamFile); PrintLowerTicks = GetNextLong(ParamFile); PrintUpperTicks = GetNextLong(ParamFile); PrintTimeLabelsXY = GetNextLong(ParamFile); PrintTimeLabelsZ = GetNextLong(ParamFile); PrintCompChar = GetNextLong(ParamFile); PrintBaseline = GetNextLong(ParamFile); PrintFieldValues = 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); FrameDistance = 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); MajorFieldPeriod = GetNextLong(ParamFile); MinorFieldPeriod = GetNextLong(ParamFile); MajorFieldLength = GetNextFloat(ParamFile); MinorFieldLength = GetNextFloat(ParamFile); FieldTickMarkWidth = GetNextFloat(ParamFile); FieldTickTop = GetNextFloat(ParamFile); FieldTickBottom = GetNextFloat(ParamFile); SkipLines(ParamFile,2); GetNextStr(ParamFile,FieldLabelFont); FieldLabelFontSize = GetNextFloat(ParamFile); FieldLabelOffset = GetNextFloat(ParamFile); SkipLines(ParamFile,2); GetNextStr(ParamFile,CompFont); CompFontSize = GetNextFloat(ParamFile); CompCharOffset = GetNextFloat(ParamFile); SkipLines(ParamFile,2); BaselineWidth = GetNextFloat(ParamFile); SkipLines(ParamFile,2); GetNextStr(ParamFile,AverageFont); AverageFontSize = GetNextFloat(ParamFile); AverageOffset = GetNextFloat(ParamFile); SkipLines(ParamFile,2); ZoomFactor = GetNextFloat(ParamFile); GridLineWidth = GetNextFloat(ParamFile); GridLineType = GetNextLong(ParamFile); GridLinePeriod = GetNextFloat(ParamFile); SkipLines(ParamFile,9); /* Read baselines if they are defined */ while (ReadLine(ParamFile,Line)) { if ((*Line != '#') && (strlen(Line) > 0)) { p = (BaselineStructPtr) malloc(sizeof(struct BaselineStruct)); strncpy(p->StationID,Line,3); p->StationID[3] = '\0'; sscanf(Line+4,"%d%d%d",&(p->Xbase),&(p->Ybase),&(p->Zbase)); p->Next = NULL; if (BaselineList == NULL) BaselineList = p; else q->Next = p; q = p; } } fclose(ParamFile); return 0; } /*--------------------------------------------------------------------------*/ /* Find if summer time is in effect in Finland at the specified time. */ /* Returns 1 if yes and 0 if not. */ /* Summer time in Finland starts on last sunday in March and ends on last */ /* sunday on October. */ /*--------------------------------------------------------------------------*/ long SummerTime(Time_sec T) { Time_struct SummerStart = {0,0,1,31, 2, 0}; /* {secs,mins,...,years} */ Time_struct SummerEnd = {0,0,1,31, 9, 0}; /* {secs,mins,...,years} */ Time_struct CurrentTime; Time_struct RefSunday = {0,0,0, 4, 0,81}; Time_sec SummerStartSec,SummerEndSec,RefSundaySec; SecsToTm(T,&CurrentTime); SummerStart.tm_year = CurrentTime.tm_year; SummerEnd.tm_year = CurrentTime.tm_year; SummerStartSec = TmToSecs(&SummerStart); SummerEndSec = TmToSecs(&SummerEnd); RefSundaySec = TmToSecs(&RefSunday); SummerStartSec -= 86400*(((SummerStartSec - RefSundaySec)/86400) % 7); SummerEndSec -= 86400*(((SummerEndSec - RefSundaySec)/86400) % 7); if ((T > SummerStartSec) && (T < SummerEndSec)) return 1; else return 0; } /*--------------------------------------------------------------------------*/ /* Create the title string using the station name, start time and end time.*/ /* If the plot contains data for a single date then the title string is of */ /* the form: NUR 1997-01-12 00-24 UT */ /* If the plot covers several days then the title string is of the form: */ /* 'NUR 1997-01-12 12 UT - 1997-01-13 12 UT' */ /*--------------------------------------------------------------------------*/ void FindTitleStr(StationPtr Station, char *Title, Time_sec StartTime, Time_sec EndTime) { char DateStr[14]; long PrintMinutes; /* Does plot start or end at full hour ? */ long OneDay; /* Does the plot fill into a single day ? */ PrintMinutes = (((StartTime % 3600) != 0) || ((EndTime % 3600) != 0)); /* First handle the station ID and start time */ strcpy(Title,Station->StationID); strcat(Title," "); if (Finnish) { /* Display time in Finnish time */ StartTime += (2 + SummerTime(StartTime))*3600; EndTime += (2 + SummerTime(StartTime))*3600; } SecsToStr(StartTime,DateStr); if (Finnish) { strncat(Title,DateStr+4,2); /* Day */ strcat(Title,"."); strncat(Title,DateStr+2,2); /* Month */ strcat(Title,"."); if (StartTime < YEAR2000) strcat(Title,"19"); else strcat(Title,"20"); strncat(Title,DateStr,2); /* Year */ } else { if (StartTime < YEAR2000) strcat(Title,"19"); else strcat(Title,"20"); strncat(Title,DateStr,2); /* Year */ strcat(Title,"-"); strncat(Title,DateStr+2,2); /* Month */ strcat(Title,"-"); strncat(Title,DateStr+4,2); /* Day */ } strcat(Title," "); strncat(Title,DateStr+6,2); /* Hour */ if (PrintMinutes) { /* add minutes */ strcat(Title,":"); strncat(Title,DateStr+8,2); } /* Handle the end time */ OneDay = ((StartTime/86400) == ((EndTime-Station->TimeStep)/86400)); SecsToStr(EndTime,DateStr); if (OneDay) { /* One day */ strcat(Title," - "); if (strncmp(DateStr+6,"0000",4) == 0) strncat(Title,"24",2); else strncat(Title,DateStr+6,2); /* Hour */ if (PrintMinutes) { /* add minutes */ strcat(Title,":"); strncat(Title,DateStr+8,2); } } else { /* Many days */ if (Finnish) { strcat(Title," - "); strncat(Title,DateStr+4,2); /* Day */ strcat(Title,"."); strncat(Title,DateStr+2,2); /* Month */ strcat(Title,"."); if (StartTime < YEAR2000) strcat(Title,"19"); else strcat(Title,"20"); strncat(Title,DateStr,2); /* Year */ } else { strcat(Title," UT - "); if (EndTime < YEAR2000) strcat(Title,"19"); else strcat(Title,"20"); strncat(Title,DateStr,2); /* Year */ strcat(Title,"-"); strncat(Title,DateStr+2,2); /* Month */ strcat(Title,"-"); strncat(Title,DateStr+4,2); /* Day */ } strcat(Title," "); strncat(Title,DateStr+6,2); /* Hour */ if (PrintMinutes) { /* add minutes */ strcat(Title,":"); strncat(Title,DateStr+8,2); } } if (!Finnish) strcat(Title," UT"); } /*--------------------------------------------------------------------------*/ /* Find the current time and put it into DateStr. */ /*--------------------------------------------------------------------------*/ void GetCurrentTime(char *DateStr) { time_t t; time(&t); strcpy(DateStr,asctime(localtime(&t))); } /*--------------------------------------------------------------------------*/ /* 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, long FrameCount) { *left = (long) (ZoomFactor*mm_to_inch72*(FrameLeft-FieldLabelOffset -6*0.35*FieldLabelFontSize)); *bottom = (long) (ZoomFactor*mm_to_inch72*(FrameBottom-FMIlogoSize-FMIlogoOffset)); *right = (long) (ZoomFactor*mm_to_inch72*(FrameLeft+FrameWidth +CompCharOffset+0.35*CompFontSize+5.0)); *top = (long) (ZoomFactor*mm_to_inch72*(FrameBottom+FrameCount*FrameHeight +(FrameCount-1)*FrameDistance+TitleOffset+0.35*TitleFontSize+5.0)); } /*--------------------------------------------------------------------------*/ /* Write the Title, Creation date and Bounding box into eps-file and do */ /* other PS-initializations. */ /*--------------------------------------------------------------------------*/ void WritePSHeader(char *StationName, char *ComponentStr) { char PSTitleStr[100]; char CreationDateStr[100]; long BBleft,BBbottom,BBright,BBtop; strcpy(PSTitleStr,StationName); strcat(PSTitleStr,"-station, "); strcat(PSTitleStr,ComponentStr); strcat(PSTitleStr,"-components"); GetCurrentTime(CreationDateStr); ComputeBoundingBox(&BBleft,&BBbottom,&BBright,&BBtop,strlen(ComponentStr)); InitializePS(PSTitleStr,"StationPlot-program",CreationDateStr, BBleft,BBbottom,BBright,BBtop); } /*--------------------------------------------------------------------------*/ /* 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; 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); } /*--------------------------------------------------------------------------*/ /* Set the MajorFieldPeriod and MinorFieldPeriod. If their values are -1 */ /* then new values will be substituted, otherwise the values read from */ /* the parameter file will be used. */ /*--------------------------------------------------------------------------*/ void SetFieldTickPeriods(long Major, long Minor) { if (MajorFieldPeriod == -1) MajorFieldPeriod = Major; if (MinorFieldPeriod == -1) MinorFieldPeriod = Minor; } /*--------------------------------------------------------------------------*/ /* Set the Major and Minor periods for the Declination axis. The values */ /* of the parameters are in arcminutes. Convert them to 2.5 arcseconds. */ /*--------------------------------------------------------------------------*/ void SetDeclTickPeriods(float MajorArcMin, float MinorArcMin, float *Major, float *Minor) { *Major = 40.0*MajorArcMin; *Minor = 40.0*MinorArcMin; } /*--------------------------------------------------------------------------*/ /* Find reasonable values for MajorFieldPeriod and MinorFieldPeriod based */ /* the total length (= Amplitude) of the vertical (= field) axis. */ /*--------------------------------------------------------------------------*/ void FindFieldTickPeriods(long Amplitude) { if (Amplitude <= 100) SetFieldTickPeriods( 50, 10); else if (Amplitude <= 400) SetFieldTickPeriods( 100, 50); else if (Amplitude <= 1000) SetFieldTickPeriods( 500, 100); else if (Amplitude <= 4000) SetFieldTickPeriods(1000, 500); else if (Amplitude <= 10000) SetFieldTickPeriods(5000,1000); else if (Amplitude <= 40000) SetFieldTickPeriods(5000,1000); else SetFieldTickPeriods(10000,5000); } /*--------------------------------------------------------------------------*/ /* Find reasonable values for MajorDeclPeriod and MinorDeclPeriod based on */ /* the total length (= Amplitude) of the vertical (= field) axis. */ /*--------------------------------------------------------------------------*/ void FindDeclTickPeriods(long Amplitude, float *Minor, float *Major) { float AmplMinute; AmplMinute = 2*Amplitude/40.0; /* Length of whole decl axis in arc minutes */ if (AmplMinute <= 5.0) SetDeclTickPeriods( 1.0, 1.0/6.0, Major, Minor); else if (AmplMinute <= 10.0) SetDeclTickPeriods( 2.0, 0.5, Major, Minor); else if (AmplMinute <= 20.0) SetDeclTickPeriods( 5.0, 1.0, Major, Minor); else if (AmplMinute <= 60.0) SetDeclTickPeriods( 10.0, 2.0, Major, Minor); else if (AmplMinute <= 120.0) SetDeclTickPeriods( 30.0, 10.0, Major, Minor); else if (AmplMinute <= 300.0) SetDeclTickPeriods( 60.0, 10.0, Major, Minor); else SetDeclTickPeriods(60.0,30.0, Major, Minor); } /*--------------------------------------------------------------------------*/ /* Draw the frame enclosing the magnetogram for a single component */ /*--------------------------------------------------------------------------*/ void DrawFrame(float yBottom) { if (PrintFrame) { LineWidthPS(FrameLineWidth); RectPS(FrameLeft,yBottom,FrameWidth,FrameHeight,0.0); } } /*--------------------------------------------------------------------------*/ /* Draw title of the plot (for current station) */ /*--------------------------------------------------------------------------*/ void DrawTitle() { float x,y; if (PrintTitle) { x = FrameLeft+FrameWidth/2; y = FrameBottom+ComponentCount*FrameHeight+(ComponentCount-1)*FrameDistance+TitleOffset; FontPS(TitleFont,TitleFontSize); TextPS(x,y,0.0,'C',TitleStr); } } /*--------------------------------------------------------------------------*/ /* Draw time labels under the time axis */ /*--------------------------------------------------------------------------*/ void DrawTimeLabels(float yBottom, Time_sec StartTime, Time_sec EndTime) { float x,dx; float TimeScale; Time_sec CurrTime; char FullDateStr[16]; char DateStr[16]; if (! PrintTimeAxis) return; FontPS(TimeLabelFont,TimeLabelFontSize); TimeScale = (FrameWidth-TimeTickLeft-TimeTickRight)/(EndTime-StartTime); x = FrameLeft+TimeTickLeft; /* First time label */ dx = 60*TimeScale*MajorTickPeriod; /* Label distance in seconds */ 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 */ if ((CurrTime==EndTime) && (strncmp(FullDateStr+6,"0000",4) == 0)) strcat(DateStr,"24"); else strncat(DateStr,FullDateStr+6,2); if (MajorTickPeriod < 60) { /* Print minutes */ strcat(DateStr,":"); strncat(DateStr,FullDateStr+8,2); } TextPS(x,yBottom-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,yBottom-TimeLabelOffset,0,'C',DateStr); } else TextPS(x,yBottom-TimeLabelOffset,0,'C',DateStr); } x += dx; CurrTime += 60*MajorTickPeriod; } while (CurrTime <= EndTime); } /*--------------------------------------------------------------------------*/ /* Draw time unit text under the last component frame */ /*--------------------------------------------------------------------------*/ void DrawTimeUnits(Time_sec StartTime, Time_sec EndTime) { float x,y; if ((! PrintTimeAxis) || (! PrintTimeLabelsZ)) return; /* 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',"Tunti (Suomen aika)"); else TextPS(x,y,0.0,'C',"Hour (UT)"); } else { /* Print day text */ if (Finnish) TextPS(x,y,0.0,'C',"Vuorokausi (Suomen aika)"); else TextPS(x,y,0.0,'C',"Day of month (UT)"); } } /*--------------------------------------------------------------------------*/ /* Draw the time axis, tick marks and vertical grid lines */ /*--------------------------------------------------------------------------*/ void DrawTimeAxis(float yBottom, 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; /* First time tick mark */ dx = 60*TimeScale*MinorTickPeriod; /* Distance between tick marks */ do { if ((CurrTime-StartTime) % (60*MajorTickPeriod) == 0) { /* Major */ if (PrintVerGridLines) { LineWidthPS(GridLineWidth); LineTypePS(GridLineType,GridLinePeriod); LinePS(x,yBottom,x,yBottom+FrameHeight); LineTypePS(0,0.0); /* Return the continuous line type */ } LineWidthPS(TickMarkWidth); if (PrintLowerTicks) LinePS(x,yBottom,x,yBottom+MajorTickLength); if (PrintUpperTicks) LinePS(x,yBottom+FrameHeight,x, yBottom+FrameHeight-MajorTickLength); } else { /* Minor */ LineWidthPS(TickMarkWidth); if (PrintLowerTicks) LinePS(x,yBottom,x,yBottom+MinorTickLength); if (PrintUpperTicks) LinePS(x,yBottom+FrameHeight,x, yBottom+FrameHeight-MinorTickLength); } x += dx; CurrTime += 60*MinorTickPeriod; } while (CurrTime <= EndTime); } /*--------------------------------------------------------------------------*/ /* Draw the field axis tick marks and horizontal grid lines and also the */ /* field values. */ /*--------------------------------------------------------------------------*/ void DrawFieldAxis(float yBottom, long BaseValue, long Amplitude) { float y,dy; float FieldScale; long CurrValue,StartValue; char FieldStr[20]; FontPS(FieldLabelFont,FieldLabelFontSize); FieldScale = (FrameHeight-FieldTickTop-FieldTickBottom)/(2*Amplitude); StartValue = BaseValue-Amplitude; CurrValue = StartValue; y = yBottom+FieldTickBottom; dy = MinorFieldPeriod*FieldScale; do { if ((CurrValue-StartValue) % MajorFieldPeriod == 0) { /* Major */ if (PrintHorGridLines) { LineWidthPS(GridLineWidth); LineTypePS(GridLineType,GridLinePeriod); LinePS(FrameLeft,y,FrameLeft+FrameWidth,y); LineTypePS(0,0.0); /* Return the continuous line type */ } LineWidthPS(FieldTickMarkWidth); LinePS(FrameLeft,y,FrameLeft+MajorFieldLength,y); LinePS(FrameLeft+FrameWidth,y, FrameLeft+FrameWidth-MajorFieldLength,y); if (PrintFieldValues) { /* Field values */ sprintf(FieldStr,"%d",CurrValue/10); TextPS(FrameLeft-FieldLabelOffset,y-0.1*FieldLabelFontSize ,0,'R',FieldStr); } } else { /* Minor */ LineWidthPS(FieldTickMarkWidth); LinePS(FrameLeft,y,FrameLeft+MinorFieldLength,y); LinePS(FrameLeft+FrameWidth,y, FrameLeft+FrameWidth-MinorFieldLength,y); } y += dy; CurrValue += MinorFieldPeriod; } while (CurrValue <= BaseValue+Amplitude); } /*--------------------------------------------------------------------------*/ /* Draw the declination axis tick marks and horizontal grid lines and also */ /* the field values. */ /*--------------------------------------------------------------------------*/ void DrawDeclinationAxis(float yBottom, long BaseValue, long Amplitude) { float y,dy; float FieldScale; float MinorDeclPeriod, MajorDeclPeriod; long CurrValue,StartValue; long Degrees, Minutes; char FieldStr[20]; FindDeclTickPeriods(Amplitude, &MinorDeclPeriod, &MajorDeclPeriod); FontPS(FieldLabelFont,FieldLabelFontSize); FieldScale = (FrameHeight-FieldTickTop-FieldTickBottom)/(2*Amplitude); /* Minor ticks */ StartValue = RoundFloat(MinorDeclPeriod*RoundFloat((BaseValue-Amplitude)/MinorDeclPeriod)); CurrValue = StartValue; y = yBottom+FrameHeight/2.0-FieldScale*(BaseValue-StartValue); dy = MinorDeclPeriod*FieldScale; do { if (CurrValue >= (BaseValue-Amplitude)) { LineWidthPS(FieldTickMarkWidth); LinePS(FrameLeft,y,FrameLeft+MinorFieldLength,y); LinePS(FrameLeft+FrameWidth,y,FrameLeft+FrameWidth-MinorFieldLength,y); } y += dy; CurrValue = RoundFloat(CurrValue+MinorDeclPeriod); } while (CurrValue <= BaseValue+Amplitude); /* Major ticks */ StartValue = RoundFloat(MajorDeclPeriod*RoundFloat((BaseValue-Amplitude)/MajorDeclPeriod)); CurrValue = StartValue; y = yBottom+FrameHeight/2.0-FieldScale*(BaseValue-StartValue); dy = MajorDeclPeriod*FieldScale; do { if (CurrValue >= (BaseValue-Amplitude)) { if (PrintHorGridLines) { LineWidthPS(GridLineWidth); LineTypePS(GridLineType,GridLinePeriod); LinePS(FrameLeft,y,FrameLeft+FrameWidth,y); LineTypePS(0,0.0); /* Return the continuous line type */ } LineWidthPS(FieldTickMarkWidth); LinePS(FrameLeft,y,FrameLeft+MajorFieldLength,y); LinePS(FrameLeft+FrameWidth,y, FrameLeft+FrameWidth-MajorFieldLength,y); if (PrintFieldValues) { /* Field values */ Degrees = CurrValue/(60*40); Minutes = (abs(CurrValue)/40) % 60; sprintf(FieldStr,"%d\\312%02d\\251",Degrees,Minutes); TextPS(FrameLeft-FieldLabelOffset,y-0.1*FieldLabelFontSize ,0,'R',FieldStr); } } y += dy; CurrValue = RoundFloat(CurrValue+MajorDeclPeriod); } while (CurrValue <= BaseValue+Amplitude); } /*--------------------------------------------------------------------------*/ /* Draw text displaying the current component to the right of the plot */ /*--------------------------------------------------------------------------*/ void DrawComponentChar(float yBottom, char Component) { float x,y; char CompStr[2] = "?"; if (!PrintCompChar) return; FontPS(CompFont,CompFontSize); x = FrameLeft+FrameWidth+CompCharOffset; y = yBottom+0.5*FrameHeight-0.1*CompFontSize; CompStr[0] = Component; TextPS(x,y,0.0,'L',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+ComponentCount*FrameHeight+(ComponentCount-1)*FrameDistance+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); } /*--------------------------------------------------------------------------*/ /* 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 != -1) Base = 10*(p->Xbase); break; case 'D' : case 'Y' : if (p->Ybase != -1) Base = 10*(p->Ybase); break; case 'Z' : if (p->Zbase != -1) Base = 10*(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) { if (!PrintBaseline) return; LineWidthPS(BaselineWidth); LineTypePS(1,8); LinePS(FrameLeft,yBase+FrameHeight/2, FrameLeft+FrameWidth,yBase+FrameHeight/2); LineTypePS(0,0); } /*--------------------------------------------------------------------------*/ /* 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,long Amplitude,Time_sec StartTime, Time_sec EndTime,long AverageTime) { float x,y,dx; float TimeScale; float FieldScale; Time_sec CurrTime; long first = 1; long DotCount = 0; long Value; LineWidthPS(DataLineWidth); LineTypePS(0,0); FieldScale = (FrameHeight-FieldTickTop-FieldTickBottom)/(2*Amplitude); TimeScale = (FrameWidth-TimeTickLeft-TimeTickRight)/(EndTime-StartTime); x = FrameLeft+TimeTickLeft; dx = AverageTime*TimeScale; CurrTime = StartTime; while (CurrTime < EndTime) { Value = ComputeAverage(Station,Component,CurrTime,AverageTime); if ((Value != MissingValue) && ((Value > BaseValue+Amplitude) || (Value < BaseValue-Amplitude))) Value = MissingValue; if (PlotMode == 1) { /* Draw continuous line between data points */ if (Value != MissingValue) { y = yBase+FieldScale*(Value-BaseValue); if (first) { MoveToPS(x,y); first = 0; } else { LineToPS(x,y); DotCount++; if (DotCount == 300) { /* This prevents stack overflow */ ClosePathPS(); /* in some 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); FillCirclePS(x,y,DotRadius,0.0,0); } } x += dx; CurrTime += AverageTime; } if (!first) ClosePathPS(); } /*--------------------------------------------------------------------------*/ /* Round the given number by specified rounding factor. */ /*--------------------------------------------------------------------------*/ long RoundToNearest(long Factor,long Number,long Direction) { if (Direction < 0) /* value < Number */ return Factor*(Number/Factor); else return Factor*(Number/Factor+1); } /*--------------------------------------------------------------------------*/ /* Get maximum value out of two or three numbers */ /*--------------------------------------------------------------------------*/ long Max2(long N1,long N2) { if (N1 == MissingValue) return 0; if (N1 > N2) return N1; else return N2; } long Max3(long N1,long N2,long N3) { long MaxValue = N1; if (N2 > MaxValue) MaxValue = N2; if (N3 > MaxValue) MaxValue = N3; return MaxValue; } /*--------------------------------------------------------------------------*/ /* Round the base line value to nearest nice value. */ /*--------------------------------------------------------------------------*/ long RoundBase(long Ave,long Amplitude,long Direction) { if (Amplitude <= 50) return RoundToNearest( 20,Ave,Direction); if (Amplitude <= 150) return RoundToNearest( 50,Ave,Direction); if (Amplitude <= 300) return RoundToNearest( 100,Ave,Direction); if (Amplitude <= 700) return RoundToNearest( 200,Ave,Direction); if (Amplitude <= 1500) return RoundToNearest( 500,Ave,Direction); if (Amplitude <= 3000) return RoundToNearest( 1000,Ave,Direction); if (Amplitude <= 7000) return RoundToNearest( 2000,Ave,Direction); if (Amplitude <= 15000) return RoundToNearest( 5000,Ave,Direction); if (Amplitude <= 30000) return RoundToNearest(10000,Ave,Direction); if (Amplitude <= 70000) return RoundToNearest(20000,Ave,Direction); if (Amplitude <=150000) return RoundToNearest(50000,Ave,Direction); return RoundToNearest(Amplitude/2,Ave,Direction); } /*--------------------------------------------------------------------------*/ /* Determine a nice value for the baseline. This is done by looking at the */ /* field amplitude and rounding the baseline value appropriately. */ /*--------------------------------------------------------------------------*/ long DetermineBase(long Max,long Min,long Ave,long MaxAmpl) { if (MaxAmpl == 0) return MissingValue; if ((Ave-Min) > (Max-Ave)) return (RoundBase(Ave,MaxAmpl,-1)); /* Round downwards */ else return (RoundBase(Ave,MaxAmpl, 1)); /* Round upwards */ } /*--------------------------------------------------------------------------*/ /* Round the total field variation amplitude upwards so that the grams */ /* will nicely fit into the figure. */ /*--------------------------------------------------------------------------*/ long RoundFieldAmplitude(long Amplitude) { if (Amplitude < 50) return 50; if (Amplitude < 100) return 100; if (Amplitude < 200) return 200; if (Amplitude < 300) return 300; if (Amplitude < 500) return 500; if (Amplitude < 1000) return 1000; if (Amplitude < 2000) return 2000; if (Amplitude < 3000) return 3000; if (Amplitude < 5000) return 5000; if (Amplitude < 10000) return 10000; if (Amplitude < 20000) return 20000; if (Amplitude < 30000) return 30000; return 50000; } /*--------------------------------------------------------------------------*/ /* Find such a field amplitude that grams for specified components will */ /* nicely fit into their boxes. */ /*--------------------------------------------------------------------------*/ long FindAmplitude(StationPtr Station,Time_sec StartTime,Time_sec EndTime,char *Components) { long Max,Min,Ampl; long Ave,Base; long MaxAmpl = 0; long i; char Comp; for (i=0;i MaxAmpl) MaxAmpl = Ampl; } return(RoundFieldAmplitude(MaxAmpl)); } /*--------------------------------------------------------------------------*/ /* Find a suitable baseline value for the field based on the maximum, */ /* minimum and average values of the field over the given period. */ /*--------------------------------------------------------------------------*/ long FindBaseline(StationPtr Station,Time_sec StartTime,Time_sec EndTime,char Component,long MaxAmpl) { long Max,Min,Ampl; long Ave; long i; FindMaxMin(Station,Component,StartTime,EndTime-StartTime,&Max,&Min); Ave = ComputeAverage(Station,Component,StartTime,EndTime-StartTime); return(DetermineBase(Max,Min,Ave,MaxAmpl)); } /*--------------------------------------------------------------------------*/ /* Draw the Finnish Meteorological Institute (FMI) logo */ /*--------------------------------------------------------------------------*/ void DrawFMILogo(void) { float xCenter; float yCenter; xCenter = FrameLeft-FMIlogoSize; yCenter = FrameBottom-FMIlogoSize/2-FMIlogoOffset; FMILogoPS(xCenter,yCenter,FMIlogoSize); FontPS(FMIlogoFont,FMIlogoFontSize); if (Finnish) TextPS(xCenter+FMIlogoSize/2+FMIlogoTextOffset, yCenter-0.1*FMIlogoFontSize, 0,'L',"Ilmatieteen laitos"); else TextPS(xCenter+FMIlogoSize/2+FMIlogoTextOffset, yCenter-0.1*FMIlogoFontSize, 0,'L',"Finnish Meteorological Institute"); } /*--------------------------------------------------------------------------*/ /* 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)); } /*--------------------------------------------------------------------------*/ /* Compute HDF from XYZ and fill the appropriate data fields in the station */ /* structure. */ /*--------------------------------------------------------------------------*/ void Compute_HDF(StationPtr s) { long *XPtr,*YPtr,*ZPtr,*HPtr,*DPtr,*FPtr; long i,count; double X,Y,Z; /* double coeff = 34377.468; = (180*60*10)/Pi */ double coeff = 137509.87; /* = 4 * (180*60*10)/Pi . Unit of D is (0.1 Arcmin)/4 = 1.5 arcseconds*/ /* Set the pointers */ XPtr = s->X; YPtr = s->Y; ZPtr = s->Z; HPtr = s->H; DPtr = s->D; FPtr = s->F; count = (s->EndTime - s->StartTime)/s->TimeStep; for (i=0;i 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); } ComponentCount = strlen(CompStr); if ((ComponentCount < 1) || (ComponentCount > 6)) { PrintUsage(argv[0],1,HeaderLineCount,UsageLineCount); return FAIL; } /*======================================*/ /* Read 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 if (!strcmp(FormatStr,"TEXT")) { status = ReadText(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 */ /*==============================================*/ /* 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; } } /*==========================================================*/ /* 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 (FieldAmplitude != 0) { FieldAmplitudeFlag = 1; FieldAmplitude *= 5; /* Convert to 0.1 nT and change definition */ /* so that field axis varies from */ /* -FieldAmplitude to +FieldAmplitude */ } /*===============================================*/ /* 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(NETWORK.StationList->StationID,CompStr); /*======================================*/ /* Go through stations and components */ /*======================================*/ Station = NETWORK.StationList; while (Station != NULL) { NewPagePS(PageNbr++); if (LandscapeMode) LandscapePS(); ZoomPS(ZoomFactor); if (strlen(HeaderStr) > 0) strcpy(TitleStr,HeaderStr); else FindTitleStr(Station,TitleStr,StartTime,EndTime); DrawTitle(); DrawAveragingValue(AverageTime); DrawTimeUnits(StartTime,EndTime); Compute_HDF(Station); if (FieldAmplitudeFlag == 0) FieldAmplitude = FindAmplitude(Station,StartTime,EndTime,CompStr); FindFieldTickPeriods(FieldAmplitude); for (CompIndex = 0; CompIndex < ComponentCount; CompIndex++) { Component = CompStr[CompIndex]; yBottom = FrameBottom+(ComponentCount-1-CompIndex)*(FrameHeight+FrameDistance); DrawFrame(yBottom); DrawComponentChar(yBottom,Component); DrawTimeAxis(yBottom,StartTime,EndTime); if ((CompIndex < ComponentCount-1) && PrintTimeLabelsXY) DrawTimeLabels(yBottom,StartTime,EndTime); if ((CompIndex == (ComponentCount-1)) && PrintTimeLabelsZ) DrawTimeLabels(yBottom,StartTime,EndTime); BaseValue = FindBaseline(Station,StartTime,EndTime,Component,FieldAmplitude); DrawBaseline(yBottom); if (Component == 'D') DrawDeclinationAxis(yBottom,BaseValue,FieldAmplitude); else DrawFieldAxis(yBottom,BaseValue,FieldAmplitude); DrawMagnetogram(Station,Component,yBottom+0.5*FrameHeight,BaseValue,FieldAmplitude,StartTime,EndTime,AverageTime); } if (FMIlogoFlag) DrawFMILogo(); ClosePagePS(); Station = Station->Next; } ClosePS(); FreeNetwork(&NETWORK); return OK; }