/**** * This files plot a series of contour lines for a 2D array of double. * 4/14/92. ****/ #include #include #include #include #include #include #include #ifndef PI #define PI 3.1415926 #endif /**** * Return true is val is between Z[i][j] and Z[i+1][j]. ****/ #define INTERSECTI(val, Z, i, j) \ (Z[i][j] <= val && val <= Z[i+1][j] || \ Z[i][j] >= val && val >= Z[i+1][j]) /**** * Return true is val is between Z[i][j] and Z[i][j+1]. ****/ #define INTERSECTJ(val, Z, i, j) \ (Z[i][j] <= val && val <= Z[i][j+1] || \ Z[i][j] >= val && val >= Z[i][j+1]) FILE *GetWriteFile(char *Ext); /* in conho.c. */ struct XY { double x, y; }; struct XYpairs { double x, y; struct XYpairs *next; }; typedef struct XYpairs *PairList; /* set up a que for printing the pairs. */ typedef PairList QDATA; struct Qlinked_list { QDATA d; struct Qlinked_list *next; }; typedef struct Qlinked_list QELEMENT; typedef QELEMENT *QLINK; typedef struct { /* a que of pairs. */ QLINK front, rear; } QUE; /* end of que definition. */ struct Isolines { double iso_val; /* z value of the contour line. */ PairList pairs; struct Isolines *next; }; typedef struct Isolines *IsoList; double ZMin(double **Z, long IXmax, long IYmax, double Dx, double Dy, struct XY * Pmin_Ptr) { double zmin; long i, j; zmin = Z[0][0]; for (i = 0; i < IXmax; i++) for (j = 0; j < IYmax; j++) if (Z[i][j] < zmin) { zmin = Z[i][j]; Pmin_Ptr->x = i * Dx; Pmin_Ptr->y = j * Dy; } return (zmin); } double ZMax(double **Z, long IXmax, long IYmax, double Dx, double Dy, struct XY * Pmax_Ptr) { double zmax; long i, j; zmax = Z[0][0]; for (i = 0; i < IXmax; i++) for (j = 0; j < IYmax; j++) if (Z[i][j] > zmax) { zmax = Z[i][j]; Pmax_Ptr->x = i * Dx; Pmax_Ptr->y = j * Dy; } return (zmax); } /* Get the isovalues from user. */ IsoList GetIsoValues(double Z_Min, double Z_Max) { IsoList head; char in_str[256]; printf("Input an isovalue or . to stop: "); scanf("%s", in_str); if (strlen(in_str) == 1 && in_str[0] == '.') return (NULL); head = (IsoList) malloc(sizeof(struct Isolines)); if (head == NULL) return (NULL); /* get the elements for the node. */ sscanf(in_str, "%lf", &head->iso_val); if (head->iso_val < Z_Min) head->iso_val = Z_Min; else if (head->iso_val > Z_Max) head->iso_val = Z_Max; head->pairs = NULL; head->next = GetIsoValues(Z_Min, Z_Max); return (head); } /**** * The isoposition is between [i][j] & [i+1][j] * * Linear interpolation is used to locate the y component of * the iso position. ****/ void IsoPositionI(PairList pair, double iso_val, double **Z, long i, long j, double Dx, double Dy) { pair->y = (j + 0.5) * Dy; if (Z[i][j] != Z[i + 1][j]) pair->x = (i + 0.5) * Dx + Dx * (iso_val - Z[i][j]) / (Z[i + 1][j] - Z[i][j]); else pair->x = (i + 1) * Dx; /* take the mid point. */ } /**** * The isoposition is between [i][j] & [i][j+1] * * Linear interpolation is used to locate the y component of * the iso position. ****/ void IsoPositionJ(PairList pair, double iso_val, double **Z, long i, long j, double Dx, double Dy) { pair->x = (i + 0.5) * Dx; if (Z[i][j + 1] != Z[i][j]) pair->y = (j + 0.5) * Dy + Dy * (iso_val - Z[i][j]) / (Z[i][j + 1] - Z[i][j]); else pair->y = (j + 1) * Dy; /* take the mid point. */ } void IsoPosition(PairList pair, double iso_val, double **Z, long i, long j, double Dx, double Dy) { if (INTERSECTI(iso_val, Z, i, j)) IsoPositionI(pair, iso_val, Z, i, j, Dx, Dy); else if (INTERSECTJ(iso_val, Z, i, j)) IsoPositionJ(pair, iso_val, Z, i, j, Dx, Dy); } PairList AllocPair(void) { PairList pair; pair = (PairList) malloc(sizeof(struct XYpairs)); if (pair == NULL) puts("...malloc error. Contour lines are not complete"); return (pair); } void GetAnIsoLine(IsoList IsoNode, double **Z, long IXmax, long IYmax, double Dx, double Dy) { long i, j; double ival = IsoNode->iso_val; PairList pair_tail; for (j = 0; j < IYmax - 1; j++) for (i = 0; i < IXmax - 2; i++) if (INTERSECTI(ival, Z, i, j) || INTERSECTJ(ival, Z, i, j)) { /* found a pair. */ if (IsoNode->pairs == NULL) { /* 1st pair. */ if (!(IsoNode->pairs = AllocPair())) return; pair_tail = IsoNode->pairs; } else { /* subsequent pairs. */ if (!(pair_tail->next = AllocPair())) return; pair_tail = pair_tail->next; } IsoPosition(pair_tail, ival, Z, i, j, Dx, Dy); } pair_tail->next = NULL; /* end of the pair list. */ } /**** * Return the quadrant. * For Frz or Arz, the peak point is usually on z axis which is * y here. We want start the isoline from the 4th quadrant. * Therefore, we use -1 for the 4th quadrant instead of 4. ****/ short GetQuadrant(double x, double y) { if (x > 0 && y >= 0) return (1); /* Include +x-axis. */ else if (x <= 0 && y > 0) return (2); /* Include +y-axis. */ else if (x < 0 && y <= 0) return (3); /* Include -x-axis. */ else if (x >= 0 && y < 0) return (-1); /* Include -y-axis. */ } /**** * Compare the angle wrt Pmax. If the angle of "This" is larger * than that of the "Next", return 1. Otherwise, return 0. * * Although atan2 returns value in the range -pi to pi, we want * to avoid it because it is slow. We compare the quadrants of * the two points first. If they are in the same quadrant, we * compare the relative positions. ****/ char OutOfOrder(PairList This, /* current pair. */ PairList Next, struct XY Pmax) { double x0, y0, x1, y1; short q0, q1; /* Quadrants. */ char out_order; if (This == NULL || Next == NULL) return (0); /* This shouldn't happen. */ x0 = This->x - Pmax.x; y0 = This->y - Pmax.y; q0 = GetQuadrant(x0, y0); x1 = Next->x - Pmax.x; y1 = Next->y - Pmax.y; q1 = GetQuadrant(x1, y1); if (q0 < q1) out_order = 0; else if (q0 > q1) out_order = 1; else { /* In the same quadrant. */ if (y0 * x1 < y1 * x0) out_order = 0; else out_order = 1; } return (out_order); } /**** * Sort the isopositions according to the angle with respect to * the position of the maximum value. ****/ void SortAnIsoLine(PairList * PairHeadPtr, struct XY Pmax) { char sorted = 0; PairList this, last; /* this=the pair being compared w/ the * next one. */ PairList sorted_head = NULL; /* the head to the sublist of * sorted pairs. */ while (!sorted && sorted_head != *PairHeadPtr) { sorted = 1; /* assume sublist is sorted. */ last = NULL; this = *PairHeadPtr; /* sublist starts at *PairHeadPtr, ends at * sorted_head. */ while (this->next != sorted_head) { /* traverse the sublist. */ if (OutOfOrder(this, this->next, Pmax)) { /* swap and move to next. */ sorted = 0; if (last == NULL) *PairHeadPtr = this->next; else last->next = this->next; last = this->next; this->next = last->next; last->next = this; } else { /* move to next. */ last = this; this = this->next; } } /* end of sublist traversing. */ sorted_head = this; } /* end of big while. */ } #if 0 /**** * Repeat the first pair of the isoline at the end of the list * so that the isoline is a loop. * Sometimes this may not be desired. ****/ void LoopAnIsoLine(PairList * PairHeadPtr) { PairList tail; tail = *PairHeadPtr; while (tail->next != NULL) tail = tail->next; tail->next = (PairList) malloc(sizeof(struct XYpairs)); if (tail->next == NULL) return; tail = tail->next; tail->x = (*PairHeadPtr)->x; tail->y = (*PairHeadPtr)->y; tail->next = NULL; } #endif void GetIsoLines(IsoList isos, double **Z, long IXmax, long IYmax, double Dx, double Dy, struct XY Pmax) { IsoList node = isos; while (node != NULL) { GetAnIsoLine(node, Z, IXmax, IYmax, Dx, Dy); SortAnIsoLine(&node->pairs, Pmax); /* LoopAnIsoLine(&node->pairs); */ node = node->next; } } char IsEmpty(QUE q) { return (q.front == NULL); } void Deque(QUE * q, QDATA * x) { QLINK temp = q->front; if (!IsEmpty(*q)) { *x = temp->d; q->front = temp->next; free(temp); } else printf("Empty que.\n"); } void Enque(QUE * q, QDATA x) { QLINK temp; temp = (QLINK) malloc(sizeof(QELEMENT)); temp->d = x; temp->next = NULL; if (IsEmpty(*q)) q->front = q->rear = temp; else { q->rear->next = temp; q->rear = temp; } } void WriteIsoLines(FILE * Isofile, IsoList IsoHead) { QUE qprint = {NULL, NULL}; /* que to be printed. */ QUE qsave = {NULL, NULL}; QDATA p; /* QDATA is PairList. */ char all_nulls; while (IsoHead != NULL) { fprintf(Isofile, "X%-8lG\tY%-8lG\t", IsoHead->iso_val, IsoHead->iso_val); Enque(&qprint, IsoHead->pairs); IsoHead = IsoHead->next; } fprintf(Isofile, "\n"); do { all_nulls = 1; while (!IsEmpty(qprint)) { Deque(&qprint, &p); if (p != NULL) { fprintf(Isofile, "%9.2lE\t%9.2lE\t", p->x, p->y); Enque(&qsave, p->next); /* add the next to a new que for a new * row. */ if (all_nulls == 1) all_nulls = 0; } else { /* fprintf(Isofile, "%9s\t%9s\t", " ", " "); */ fprintf(Isofile, "\t\t"); Enque(&qsave, p); /* add a null to the new que. */ } } fprintf(Isofile, "\n"); qprint = qsave; qsave.front = NULL; qsave.rear = NULL; } while (!all_nulls); } void FreePairs(PairList Pairs) { PairList this_pair; while (Pairs != NULL) { this_pair = Pairs; Pairs = Pairs->next; free(this_pair); } } void FreeIsoLines(IsoList IsoHead) { IsoList this_iso; PairList pair_head; while (IsoHead != NULL) { FreePairs(IsoHead->pairs); this_iso = IsoHead; IsoHead = IsoHead->next; free(this_iso); } } void IsoPlot(double **Z, /* the 2D array Z[i][j]. */ long int IXmax, long int IYmax, /* the 0<=i<=IXmax, 0<=j<=IYmax. */ double Dx, double Dy) { /* the gridline separations. */ FILE *isofile; /* send isolines to this file. */ struct XY pmin, pmax; /* the xy positions of the min & max. */ double zmin, zmax; IsoList isos; char fname[128] = "iso"; /* Locate the min & max. */ zmin = ZMin(Z, IXmax, IYmax, Dx, Dy, &pmin); zmax = ZMax(Z, IXmax, IYmax, Dx, Dy, &pmax); if (zmin == zmax || IXmax * IYmax < 4) { printf("...Not enough data for contour plot\n"); return; } isofile = GetWriteFile(fname); if (isofile == NULL) return; printf("The range of the value is %lf to %lf.\n", zmin, zmax); isos = GetIsoValues(zmin, zmax); GetIsoLines(isos, Z, IXmax, IYmax, Dx, Dy, pmax); WriteIsoLines(isofile, isos); FreeIsoLines(isos); fclose(isofile); }